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ABSTRACT 


Hybrid  state  estimation  problems  are  statistical  estimation  problems  in 
which  both  continuous-valued  and  discrete-valued  states  and  parameters  occur. 
The  hybrid  state  model  provides  both  a  natural  formulation  for  many  types  of 
surveillance  and  tracking  problems  and  a  powerful  framework  for  deriving  theo 
retically  optimal  and  practical  suboptimal  tracking  algorithms.  This  report 
describes  research  in  developing  tools  for  evaluating  the  performance  of  opti 
mal  hybrid  state  estimation  problems. 


s.  -r.  f 
r 


ALPHATECH,  INC 


ACKNOWLEDGEMENT 


This  research  was  supported  by  Code  340R  of  the  Naval  Air  Systems  Command 
under  contract  N00014-84-C-03078  with  the  Office  of  Naval  Research. 


P^esfon^ToT 
hj>s  CRA&T 
I  nT‘C  TAB 

ililir 

i  By . 

ibutior,  I 


Dht  /  *’d/oT 


ALPHATECH,  INC 


ALPHATECH,  INC 


CONTENTS  (Continued) 


MONTE-CARLO  PERFORMANCE  ANALYSIS  BASED  ON  THE  REPRESENTATION 
THEOREM . 

6.1  INTRODUCTION . 

6.2  GALDOS  RATE  DISTORTION  METHOD  . 

6.3  DIRECT  MONTE-CARLO  PERFORMANCE  ANALYSIS  . 

6.4  CONCLUSIONS  . 

A  RANDOM  POINT  PROCESS  APPROACH  TO  MULTIOBJECT  TRACKING.  .  . 


7.1  INTRODUCTION . 

7.2  RANDOM  POINT  PROCESS  MODELS  .  .  .  . 

7.3  LIKELIHOOD  RATIO  FUNCTIONS . 

7.4  ESTIMATION  ALGORITHM . 

7.5  APPLICATION  TO  PERFORMANCE  ANALYSIS 

7.6  CONCLUDING  REMARKS . 


NUMERICAL  EXAMPLES  AND  CONCLUSIONS 


8.1  NUMERICAL  COMPARISON  OF  PERFORMANCE  METHODS  . 

8.2  CONCLUDING  REMARKS . 

REFERENCES  . 

APPENDIX  A  -  SECOND-ORDER  EXPANSION  OF  MINIMUM  MEAN  SQUARE  ERROR  .  .  . 


ALPHATECH,  INC 


SECTION  1 
INTRODUCTION 

Although  the  mathematical  structure  of  the  optimal  solution  of  the  multi¬ 
object  tracking  problem  is  well  understood,  the  excessive  computational  com¬ 
plexity  of  the  optimal  tracking  algorithm  prohibits  its  practical  realization 
using  even  the  largest  and  fastest  computers  available  in  the  foreseeable 
future.  For  this  reason,  many  different  tracking  algorithms  have  been  devel¬ 
oped  which  sacrifice  optimal  performance  for  the  sake  of  computational  feasi¬ 
bility.  Intuition  suggests  that  some  suboptimal  approximations  to  the  optimal 
algorithm  result  in  little  performance  loss,  but  there  exists  no  reliable 
quantitative  analysis  of  suboptimal  performance  that  shows  this  is,  in  fact, 
the  case. 

To  a  limited  extent,  it  is  possible  to  compare  one  suboptimal  algorithm 
to  another  suboptimal  algorithm  by  numerically  simulating  the  performance  of 
both.  Such  simulations  provide  useful  information  about  the  relative  perfor¬ 
mance  of  different  suboptimal  algorithms,  but  they  do  not  give  reliable  quan¬ 
titative  measures  of  performance  without  extensive  simulation.  In  all  but  the 
simplest  multiobject  tracking  scenarios,  it  is  expensive  and  time-consuming 
to  compute  an  adequate  number  of  simulation  samples  for  reliable  Monte-Carlo 
analysis.  For  this  reason,  it  is  difficult  to  use  simulation  to  investigate 
how  an  algorithm's  performance  depends  on  algorithmic  and  environmental 
parameters. 
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While  simulating  relative  performance  of  suboptimal  algorithms  is  diffi¬ 
cult,  simulating  optimal  performance  to  determine  the  absolute  performance  of 
a  suboptimal  algorithm  is  effectively  impossible.  In  all  but  the  simplest 
cases,  it  is  computationally  impossible  to  compute  even  one  simulation  sample 
of  the  optimal  tracking  algorithm,  let  alone  an  adequate  number  of  samples  for 
Monte-Carlo  analysis.  Although  one  can  obtain  some  insight  by  comparing  simu¬ 
lated  performance  of  optimal  and  suboptimal  algorithms  for  a  small  number  of 
very  limited  cases,  it  is  not  possible  to  obtain  quantitative  measures  of  per¬ 
formance  in  this  way. 

What  are  lacking  are  analytical  methods  for  analyzing  the  quantitative 
performance  of  optimal  and  suboptimal  multiobject  tracking  algorithms.  Because 
of  the  excessive  computational  requirements  of  simulating  such  algorithms, 
there  is  no  feasible  alternative  to  analytical  methods  of  performance  analysis. 
Such  analytical  methods  do  exist  and  have  proven  useful  for  analyzing  per¬ 
formance  in  other,  much  simpler  types  of  hypothesis  testing  and  estimation 
problems.  There,  efficient  analytical  performance  analysis  techniques  have 
greatly  facilitated  the  design  of  practical  algorithms.  The  basic  research 
problem  here  is  to  find  similarly  effective  analytical  techniques  for  the  more 
difficult  hypothesis  testing  and  estimation  problems  of  multiobject  tracking. 

We  base  our  approach  to  this  problem  on  studying  multiobject  tracking 
problems  as  a  special  case  of  the  general  class  of  hybrid  state  estimation 
problems  [1]  , [2]*.  Hybrid  state  estimation  problems  are  statistical  estima¬ 
tion  problems  in  which  both  continuous-valued  and  discrete-valued  states  and 
parameters  occur.  Thus,  hybrid  state  estimation  includes  both  classical 

*Reference  are  indicated  by  numbers  in  square  brackets.  The  list  appears  at 
the  end  of  the  main  body  of  this  report. 
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parameter  estimation  and  hypothesis  testing  as  special  cases.  In  our  previous 
work  [1],[2],  we  have  found  that  the  hybrid  state  model  provides  both  a  nat¬ 
ural  formulation  for  all  types  of  multiobject  tracking  problems  and  a  powerful 
framework  for  deriving  optimal  and  suboptimal  tracking  algorithms.  The  im¬ 
portance  of  the  hybrid  state  system  viewpoint  in  the  research  described  here 
is  that  it  gives  us  a  bridge  between  the  simpler  estimation  and  hypothesis 
testing  problems  for  which  effective  performance  analysis  techniques  exist 
and  the  more  difficult  estimation  and  hypothesis  testing  problems  of  multi¬ 
object  tracking  for  which  effective  techniques  have  yet  to  be  developed. 

Our  overall  approach  is  to  consider  a  hierarchy  of  hybrid  state  estima¬ 
tion  problems  of  increasing  difficulty,  starting  at  relatively  simple  problems 
whose  solutions  are  readily  obtained  and  leading  to  more  difficult  multiobject 
tracking  problems.  To  carry  this  out,  it  is  convenient  to  study  three  dif¬ 
ferent  types  of  hybrid  state  estimation  problems.  In  the  first  type  (Type  I) 
of  problem  we  are  interested  only  in  estimating  continuous-valued  states  and 
discrete-valued  states  enter  the  problem  only  as  measurement  noise.  Type  I 
problems  include  those  in  which  the  number  of  targets  is  known  but  the  associ¬ 
ation  of  measurements  with  targets  and  false  alarms  is  not  known.  Since  Type 
I  problems  are  continuous  state  estimation  problems,  one  can  try  to  extend  many 
different  available  methods  of  performance  analysis  to  this  type  of  problem. 
What  is  unconventional  about  such  problems  (and  tends  to  make  performance 
analysis  difficult)  is  that  the  discrete-valued  noise  is  non-Gaussian  (being 
a  point  process  in  nature)  and  nonadditive  (often  being  multiplicative). 

The  second  type  (Type  II)  of  hybrid  state  estimation  problem  generalizes 


the  Type  I  problem  in  the  following  way.  In  this  problem  we  are  still  inter¬ 
ested  in  estimating  continuous  states,  but  the  discrete  states  enter  into  the 
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problem  In  a  much  more  fundamental  way  than  In  Type  1  problems.  In  partic¬ 
ular,  the  discrete  states  are  assumed  to  form  a  Markov  process.  Note  that  we 
assume  that  the  continuous  states  cannot  affect  the  discrete  states,  although 
the  discrete  states  can  affect  the  continuous  states.  Problems  of  tracking 
targets  which  can  maneuver,  or  estimating  parameters  in  a  system  whose  dy¬ 
namics  change  abruptly  (e.g.,  due  to  failures)  can  be  formulated  as  Type  II 
hybrid  state  estimation  problems.  Finally,  we  identify  a  third  type  of  hybrid 
state  estimation  problem  (Type  III)  which  are  the  same  as  Type  II  problems 
except  that  the  estimation  problem  includes  the  estimation  of  discrete  states 
as  well  as  continuous  states.  Type  III  problems  include  the  most  realistic 
and  most  difficult  multiobject  tracking  problems  in  which  the  number  of  tar¬ 
gets  is  unknown  and  one  is  interested  in  both  tracking  and  detecting  targets 
entering  and  leaving  a  surveillance  area. 

This  report  describes  the  results  of  our  research  in  analyzing  the  per¬ 
formance  of  hybrid  state  estimation  problems.  Section  2  describes  the  gen¬ 
eral  hybrid  state  model  of  interest  and  several  concrete  examples  relevant 
to  tracking  problems;  It  also  defines  a  classification  (Types  I,  II,  III)  of 
hybrid  state  estimation  problems.  Section  3  analyzes  a  simple  hybrid  state 
estimation  problem  in  detail  -  indicating  the  difficulties  inherent  in  deter¬ 
mining  optimal  performance  for  this  class  of  problems.  Section  4  describes  a 
Cramer-Rao  lower  bound  on  mean  square  error  for  a  large  class  of  Type  I  hybrid 
state  problems.  In  Section  5  we  discuss  a  rate  distortion  approach  for  Type 
II  problems,  and  in  Section  6  we  describe  a  Monte-Carlo  approach  based  on  the 
representation  theorem  which  applies  to  very  general  hybrid  state  problems. 

The  results  of  Sections  3  through  6  apply  to  finite  dimensional  hybrid  state 
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problems  described  in  Section  2.  Extension  to  hybrid  state  problems  of  multi¬ 
object  tracking  is  the  subject  of  Section  7.  There  we  introduce  a  random 
point  process  model  of  multiobject  tracking  and  discuss  its  application  to 
tracking  and  performance  analysis.  Section  8  concludes  the  report  and  dis¬ 
cusses  some  further  research  problems. 
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SECTION  2 

HYBRID  STATE  ESTIMATION 

2.1  GENERAL  HYBRID  STATE  MODEL 

Multiobject  tracking  problems  involve  both  the  estimation  of  continuous¬ 
valued  parameters  such  as  target  positions  and  velocities,  and  the  the  testing 
of  discrete  hypotheses  such  as  the  association  of  measurement  data  with  targets. 
Thus,  it  is  natural  to  pose  the  problem  of  multiobject  tracking  as  a  hybrid 
state  estimation  problem  —  that  is,  estimation  for  a  partially  observed  Markov 
process  with  discrete-  and  continuous-valued  states.  In  our  previous  work  [1], 
[2]  we  have  shown  that  the  hybrid  state  viewpoint  provides  a  powerful  frame¬ 
work  in  which  to  formulate  all  types  of  multiobject  tracking  problems.  The 
research  reported  here  is  continuing  this  approach  by  developing  performance 
analysis  techniques  for  several  types  of  hybrid  state  estimation  problems 
relevant  to  multiobject  tracking. 

To  focus  our  research,  we  will  consider  a  particular  class  of  hybrid 
state  system,  which  we  call  Gauss-Markov  hybrid  state  systems.  This  is  the 
class  of  systems  which  can  be  described  by  the  equations 

x(t+l)  »  A(q( t)  )x( t)  +  B (q( t) )w( t)  (2-la) 

y(t)  -  C (q ( t )  )x( t )  +  D(q(t))v(t)  (2-lb) 

where  the  variables  are  defined  so  that 
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q{t)  =  finite  state  Markov  chain  with  time  invariant 
transition  probabilities; 

x(t)  =  n-dimensional  state  process; 

w(t)  =  p-dimensional  white  Gaussian  state  process  noise; 

y(t)  =■  m-dimensional  measurement  process; 

v(t)  =*  r-dimensional  white  Gaussian  measurement  noise. 

For  a  fixed  value  of  q  of  the  Markov  chain,  A(q),  B(q),  C(q) ,  and  D(q)  are 
matrices  with  dimensions  corresponding  to  the  vectors  described  above.  One 
important  point  to  note  about  the  model  described  by  Eq.  2-1  is  that  the  dis¬ 
crete  state  process  q(t)  may  drive  the  continuous  state  x(t)  but  not  conversely 
This  assumption  simplifies  analysis  and  still  permits  the  treatment  of  many 
(but  not  all,  see  [2])  realistic  multiobject  tracking  problems. 

We  have  assumed  that  the  state  and  measurement  equations  are  linear  in 
the  continuous  state  variable  x  in  order  to  focus  on  the  essential  hybrid 
state  aspects  of  our  problem.  Some  of  the  methods  we  discuss  will  extend  to 
nonlinear  hybrid  state  problems  of  the  general  form 

x(t+l)  =  f (x(t) ,q(t) )  +  B(q(t))w(t)  (2-2a) 

y(t)  =  h(x(t) ,q(t) )  +  D(q(t))v(t)  (2-2b) 

where  B(q)  and  D(q)  are  matrices  as  above,  but  f(x,q)  and  h(x,q)  are  now 
allowed  to  be  nonlinear  functions  of  x  for  each  fixed  value  of  q.  The  non¬ 
linear  model  of  Eq.  2-2  allows  the  treatment  of  range-azimuth  measurement 
nonlinearities  and  passive  tracking  measurement  nonlinearities,  but  the  re¬ 
sulting  performance  analysis  problem  will  include  all  the  usual  difficulties 
of  analyzing  nonlinear  estimation  problems  in  addition  to  the  difficulties  of 
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analyzing  hybrid  state  problems.  In  the  next  subsection  we  consider  several 
examples  of  Gauss-Markov  hybrid  state  models  relevant  to  tracking. 


2.2  TRACKING  EXAMPLES  OF  HYBRID  STATE  MODELS 

We  will  give  four  examples  in  this  subsection  to  illustrate  how  one  may 
use  the  Gauss-Markov  hybrid  state  equations  of  2-1  to  model  tracking  problems. 
These  models  are  simplifications  of  real  tracking  problems;  in  particular, 
they  approximate  complex  transducers  and  signal  processors  with  simple  hybrid 
state  statistical  models.  We  have  chosen  this  type  of  approximation  to  focus 
our  attention  on  the  "backend”  tracking  problem  (which  is  the  object  of  our 
research)  and  to  avoid  less  important  details  of  “frontend"  processing. 

The  first  example  is  a  simple  model  of  a  sensor  and  a  single  target.  The 
sensor  produces  one  measurement  per  time  period;  each  measurement  is  either  a 
measurement  of  the  target's  state  together  with  background  noise,  a  measure¬ 
ment  of  independent  random  clutter,  or  a  measurement  of  ambient  background 
noise  alone.  The  target  moves  according  to  a  simple  linear  Gaussian  stochas¬ 
tic  difference  equation.  The  hybrid  state  model  for  this  example  is  given  as 


follows. 


xi(t+l)  =  xi(t)  +  X2(t) 


X2(t+1)  *>  X2(t)  +  w(t) 


y(t)  -  qi(t)[xi(t)  +  Vi(t)]  +  (l-qi(t) )q2(t)v2(t) 


(2-3a) 


(2-3b) 


(2-4) 


Eq.  2-3  is  the  state  equation  and  is  the  usual  random  acceleration  model  in 
which  xi  denotes  the  target  position  and  X2  denotes  the  target  velocity.  The 
white  noise  w  is  thus  an  assumed  random  acceleration.  Eq.  2-4  is  the  measure¬ 
ment  equation  for  this  example  and  contains  all  the  discrete  state  dependences 
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We  assume  that  qi(t)  and  q2<t)  are  independent  discrete  state  random  variables 
taking  variables  taking  the  values  0  or  1.  Assume  that  qi(t)  equals  1  with 
probability  pj.  Then  for  each  time  period  t,  the  measurement  y(t)  is  either 

a  target  measurement  x^(t)  with  some  error  vi(t)  with  probability  pj_; 
a  clutter  measurement  V2(t)  with  probability  (l-pi)P2; 
or  the  0  measurement  with  probability  (l~Pi)(l~P2) • 

Thus,  pi  is  the  probability  of  detecting  the  target,  ( l~Pi)P2  is  the  proba¬ 
bility  of  detecting  a  false  alarm,  and  (l-pi)(l-p2)  is  the  probability  of 
detecting  nothing  at  all.  These  probabilities  might  be  selected  to  correspond 
to  the  receiver  operating  characteristics  of  a  particular  sensor  and  signal 
processor  frontend.  The  problem  of  this  example  is  to  estimate  xi(t)  and 
X2(t),  the  target  position  and  velocity  states.  As  we  will  see  below,  the 
discrete  variables  qi(t)  and  q2(t)  are  not  desired  and  do  not  have  to  be  esti¬ 
mated  in  this  example  although  a  practical  tracking  algorithm  might  in  fact 
have  a  detector  that  would  estimate  these  discrete  variables.  This  point  will 
be  discussed  in  the  next  subsection  where  we  classify  hybrid  state  estimation 
problems  more  precisely  into  three  types  (Types  I,  II,  and  III). 

The  second  example  is  a  variation  of  the  first  to  include  data  associa¬ 
tion  ambiguities  in  the  measurement  model.  In  this  example  there  are  two 
targets  with  states  xi,x2  and  X3,X4,  respectively.  The  sensor  observes  two 
returns  yi  and  y2,  but  it  does  not  observe  the  origin  of  the  measurements. 

The  Gauss-Markov  hybrid  state  model  is  given  by 

xi(t+l)  -  Xi(t)  +  x2(t)  (2-5a) 

x2( t+1)  -  X2(t)  +  wj^t)  ( 2-5b) 
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x3(t+l)  *  X3(t)  +  X4(t)  ( 2-5c) 

x4( t+1)  *  X4(t)  +  W2(t)  .  (2-5d) 

Yl(t)  =  q(t)x1(t)  +  (l-q(t) )x2(t)  +  v3(t)  (2-6a) 

Y2(t)  =  (l-q(t)  )x3(t)  +  q(t)x2(t)  +  v2(t)  (2-6b) 


Equation  2-5  is  the  usual  state  model  for  two  independent  targets.  In 
Eq.  2-6  we  assume  that  q(t)  takes  the  values  0  or  1,  both  with  probability 
1/2.  Thus,  the  sensor  measures  both  target  positions,  but  the  ordering  of 
these  measurements  is  randomly  permuted. 

The  third  example  is  an  extension  of  the  first  to  allow  the  target  to 
make  significant  maneuvers  at  unpredictable  times.  The  measurement  model  Eq. 
2-4  remains  the  same  ,but  the  state  model  becomes  the  hybrid  state  equation 

xi(t+l)  -  x3(t)  +  x2(t)  (2-7a) 

x2(t+l)  -  x2(t)  +  q3(t)w(t)  (2-7b) 

In  Eq.  2-7b  the  discrete  state  process  q3(t)  is  a  Markov  chain  taking  two 
values  Qi  and  Q2  with  transition  probabilities  P^j  of  jumping  from  discrete 
state  Qi  to  state  Qj.  The  new  discrete  state  models  two  types  of  random 
acceleration  distinguished  by  different  covariances.  For  example,  the  covar¬ 
iance  Qi  might  have  a  small  value  representing  small  perturbations  of  essen¬ 
tially  constant  velocity  motion  and  the  covariance  Q2  might  have  a  large 
value  representing  occasional  significant  target  maneuvers.  By  appropriately 
selecting  the  transition  probabilities  Pjj,  one  can  adjust  the  average  time 
between  maneuvers.  Note  that  the  problem  in  this  example  is  to  estimate  x3(t) 
and  x2(t).  It  may  be  necessary  to  estimate  q3(t)  also,  but  that  is  only 
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incidental  to  the  main  tracking  objective.  In  a  sense,  qj(t)  is  a  nuisance 
parameter  as  far  as  tracking  is  concerned.  Note  that  in  a  different  version 
of  this  problem,  q3(t)  might  also  be  desired  (e.g.,  this  would  be  the  case  if 
detecting  maneuvers  implied  significant  tactical  information  such  as  that  a 
missile  launch  was  about  to  occur) . 

The  final  example  is  another  extension  of  the  first  example  in  which  we 
now  model  the  appearance  and  disappearance  of  the  single  target.  Let  q3(t) 
denote  a  Markov  chain  taking  two  values  0  and  1  with  transition  probabilities 
Pij  of  jumping  from  discrete  state  1  to  state  j.  When  q3(t)  ■  1,  we  consider 
that  the  target  is  present  and  when  q3(t)  -  0,  we  consider  that  the  target  is 
absent.  The  state  equation  for  this  example  is  the  same  as  Gq.  2-3;  the  mea¬ 
surement  equation  becomes  the  following. 


y(t)  -  qi(t) [q3(t)xi(t)  +  v^(t)  ]  +  (l-qi(t) )q2(t)v2(t) 


(2-8) 


The  problem  in  this  example,  unlike  the  preceding  two  examples,  is  to  estimate 
q3(t),  as  well  as  x^(t)  and  x2(t).  That  is,  it  is  desired  to  know  whether 
the  target  is  present  (i.e.,  observable)  as  well  as  what  are  its  position  and 
velocity.  Indeed,  in  some  cases  the  position  and  velocity  may  only  be  rele¬ 
vant  if  q3(t)  ■  1. 

2.3  CLASSIFICATION  OF  HYBRID  STATE  ESTIMATION  PROBLEMS 

It  is  convenient  to  define  three  types  of  hybrid  state  estimation  prob¬ 
lems  based  partly  on  the  structure  of  the  hybrid  state  model  and  partly  on  the 
desired  estimation  criterion.  Type  I  problems  are  modeled  by  hybrid  state 
equations  of  the  form  of  Eq.  2-1  where  the  discrete  state  process  q(t)  is  a 
sequence  of  independent  finite  state  random  variables.  The  estimation  problem 
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is  Co  estimate  the  state  x(t)  given  measurements  y(s)  up  to  time  t.  The 
important  point  to  note  is  that  q(t)  plays  the  role  of  a  nuisance  parameter 
or  an  additional  noise  process  analogous  to  w(t)  and  v(t).  It  is  not  neces¬ 
sary  to  estimate  q(t)  any  more  than  it  is  necessary  to  estimate  w(t)  or  v(t) 
to  obtain  an  estimate  of  x(t).  This  is  not  to  say  that  some  estimation  algo¬ 
rithm  may,  in  fact,  estimate  q(t)  as  an  intermediate  step  in  obtaining  an 
estimate  of  x(t). 

The  significance  of  Type  I  problem  lies  in  the  fact  that  they  are  essen¬ 
tially  continuous  state  estimation  problems  and  therefore,  many  classical  per 
formance  analysis  techniques  (such  as  Cramer-Rao  type  methods)  apply  to  this 
type  of  problem,  at  least  in  theory.  As  we  will  see  later  in  this  report, 
there  is  a  considerable  difficulty  in  applying  the  theory  because  the  distri¬ 
butions  are  non-Gaussian  (being  a  mixture  of  Gaussian  distributions). 

We  will  often  consider  a  subclass  of  Type  I  problems  for  which  the  con¬ 
tinuous  state  dynamics  are  independent  of  the  discrete  variables,  namely 
Eq.  2-1  is  replaced  by 


x(t+l)  -  Ax( t)  +  Bw(t) 


y(t)  -  C(q(t)  )x(t)  +  D[q(t))v(t) 


(2-9a 


(2-9b 


Even  within  this  subclass  of  Type  I  problems  it  is  possible  to  consider  many 
different  tracking  examples  of  interest.  The  first  and  second  examples  of 
the  previous  subsection  are  members  of  this  class  of  Type  I  problem.  One  can 
easily  generalize  these  examples  to  model  problems  with  multiple  targets  and 
multiple  returns  per  time  period,  provided  that  the  number  of  targets  is  con¬ 
stant  and  known  (i.e.,  detecting  a  target's  presence  or  absence  is  not  a  part 
of  the  problem) . 


ALPHATECH,  INC 


Type  II  problems  are  given  by  the  same  equations  (namely  Eq.  2-1),  but 
the  discrete  state  process  q(t)  is  allowed  to  be  a  finite  state  Markov  chain 
rather  than  a  sequence  of  independent  finite  state  random  variables.  Never¬ 
theless,  it  is  still  desired  to  estimate  only  the  continuous  state  x(t)  in 
Type  II  problems,  although  even  more  so  than  in  Type  I  problems  will  it  be 
necessary  to  estimate  the  discrete  state  q(t)  as  an  intermediate  step  in 
obtaining  the  estimate  of  x(t).  The  second  example  of  tracking  a  maneuvering 
target  is  a  Type  II  hybrid  state  estimation  problem. 

Let  us  make  the  distinction  between  Type  I  and  II  clear.  In  Type  I  prob¬ 
lems  the  joint  process  <x(t),y(t)>  is  Markov.  In  Type  II  problems  this  is 
no  longer  true;  the  Markov  process  is  the  joint  process  <x(t),y(t),q(t)>  (or 
at  least  some  component  of  q(t)).  This  property  of  Type  II  problems  compli¬ 
cates  performance  analysis  considerably.  However,  the  estimation  criterion 
depends  just  on  x(t);  for  example,  it  is  a  mean  square  error  criterion  for 
estimating  x(t).  This  assumption  at  least  gives  us  a  simple  measure  of  per¬ 
formance  to  consider. 

A  subclass  of  Type  II  problems  has  the  property  that  q(t)  is  equal  to 
q(0)  for  all  times  t.  These  are  the  so-called  multiple  model  problems  some¬ 
times  used  in  adaptive  estimation.  We  will  not  study  this  class  of  problems 
and  refer  to  reader  to  [3]  for  further  information. 

Type  III  problems  generalize  the  Type  II  problem  by  allowing  the  estima¬ 
tion  criterion  to  depend  on  both  x(t)  and  on  q(t).  This  type  of  hybrid  state 
estimation  problem  also  has  the  model  of  Eq.  2-1  in  which  q(t)  is  allowed  to 
be  a  general  finite  state  Markov  chain.  The  third  example  of  tracking  a  single 
target  which  appears  and  disappears  at  random  times  is  a  Type  III  problem  if 
we  are  interested  in  knowing  when  the  target  is  present  and  when  it  is  absent. 
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Type  III  problems  include  the  most  realistic  multiobject  tracking  prob¬ 
lems  in  which  the  total  number  of  targets  is  unknown  and  in  which  one  is 
interested  in  some  aspect  of  the  identity  of  the  target  (if  only  to  know  that 
an  object  being  tracked  at  the  current  time  is  the  same  as  an  object  that  was 
being  tracked  at  an  earlier  time;  that  is,  if  track  continuity  is  an  issue). 
In  such  problems  it  is  difficult  to  define  meaningful  measures  of  performance 
especially  if  they  are  to  be  amenable  to  analysis.  In  the  rest  of  this  paper 
we  will  discuss  performance  analysis  for  different  hybrid  state  estimation 
problems . 
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SECTION  3 


ANALYSIS  OF  A  SIMPLE  PROBLEM 


3.1  INTRODUCTION 

A  simple  hybrid  state  estimation  problem  of  interest  is  described  by  the 


measurement  equation 


y  =  x  +  qvi  +  (l-q)vo 


(3-D 


where  x,  vj_  and  vg  are  normal  random  variables  and  q  -  0  or  1.  Assume  that 
x,  vl>  v0  and  <1  are  independent  and  assume  tht  the  distributions  are 

x  ~  N(x,  o2) 

vj  ~  N  (O,  0^2  j 

vo  ~  N(0,  oq2) 

P{q=l}  *  e 


We  are  interested  in  approximating  the  minimum  mean  square  error  of  estimating 
x  given  the  measurement  y.  Equation  3-1  presents  an  example  of  estimation 
with  non-Gaussian  measurement  error.  This  is  true  in  general  for  hybrid  state 
estimation  problems. 


The  measurement  error 


qvi  +  (l-q)vQ 


is  a  random  mixture  of  two  different  random  variales  (v^  and  vo).  For  example 
this  model  might  represent  the  problem  of  estimating  x  given  a  probability  e 
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of  taking  a  bad  measurement  (with  large  variance  c»i2)  and  probability  1  -  e  of 
taking  a  good  measurement  (with  small  variance  oo^). 

Suppose  x  is  the  conditional  mean  of  x  given  by 

x  =■  E  (x  |  y  }  . 

Then  the  minimum  mean  square  error  V  is  given  by 

V  =  E{(x-x)2}  . 

Our  general  goal  is  to  develop  good  methods  for  approximating  V  in  hybrid 
state  estimation  problems.  In  this  section  we  discuss  several  approximations 
to  develop  some  feeling  for  different  methods'  computability  and  accuracy. 

In  the  following  subsections  we  consider  direct  calculation  of  V,  asymptotic 
expansion  of  V  (for  small  e) ,  the  Cramer-Rao  lower  bound  of  V,  and  a  rate 
distortion  lower  bound  of  V.  We  also  consider  numerical  examples  of  the  dif¬ 
ferent  approximations  to  study  their  relative  accuracy. 

3.2  ASYMPTOTIC  EXPANSION  OF  OPTIMAL  PERFORMANCE 
3.2.1  Exact  Expression  of  Minimum  Mean  Square  Error 

For  the  class  of  Gauss-Markov  hybrid  state  systems  defined  in  Section  2, 
it  is  easy  to  describe  the  minimum  mean  square  estimator;  namely  the  condi¬ 
tional  mean.  For  this  simple  problem  we  have 


x(y)  -  E (x | y,  q=0 }  P{q=Ojy}  +  E{x|y,  q-1 }  P{q-l|y} 

where  E{x|y,  q=0}  and  E{x|y,  q=l }  are  the  solution  of  linear  least  squares 
problems,  and  Pjq*0|y}  and  p{q=l|y}  are  easily  obtained  from  Bayes'  rule.  For 
the  problem  at  hand 
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E {x | y »  q-i} 


o l*-  x  +  o*y 
+  a2 


[2 x(o2  +  Oi2) f  1^2  exp  {  -  1/2  (y~x)2(o2  +  Oi2)'1  }  P{q-i} 


-  1/2  _ 

[2x(o2  +  012)  ]  exp  {  -  1/2  (y-x)2(o2  +  °02)-M  p{q=1} 

-  1/2  _ 

+  [2i Ka2  +  oq2)  ]  exp {  -  1/2  (y-x)2 

(a2  +  oq2)-1}  P{q=*0} 


Thus,  we  can  write  x  as 


X  ■  X  +  5 


(l-e)bo  $o(D  +  e  bi 

(1-0  <|>o(0  +  e  *l(0 


(3-2) 


where 


y-x 


bi  ■  o2(o2  +  Oi2) 


(3-3) 


4>i(0  -  [2ir(  o2  +  oi2)]"  1/2  exp  {  -  1/2  S2(o2  +  Oi2)-* }  (3-4) 

However,  we  are  not  interested  in  x  itself,  but  in  its  mean  square  performance 


V.  This  is  found  from 


V  ■  E  {(x-x)2  }  -  E  {(x-x)2  } 
■  o2  -  I 


(3-5) 


where  I  =*  e{(x-x)2}  is  the  integral 


•  [(l-e)bo  4>o(  ^)  +  e  bl  $l(0]2 

I  -  /  52  - 

1-  (l-e)taU)  +  e  <fri(0 


(3-6) 
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Even  for  this  simple  problem,  the  integral  I  does  not  have  a  closed  form  solu¬ 
tion  and  we  must  resort  to  numerical  or  analytic  approxmation  of  it.  This 
type  of  integral,  in  which  a  sum  of  Gaussian  distribution  appears  as  a  denom¬ 
inator  of  the  integrand,  is  common  to  all  Gauss-Markov  hybrid  state  problems. 
It  occurs  in  the  exact  expression  of  minimum  mean  square  error  and  in  the  ex¬ 
pression  for  the  Cramer-Rao  bound  as  we  will  show  later.  Thus,  it  is  impor¬ 
tant  that  we  learn  to  approximate  such  integrals  well.  In  the  next  subsection 
we  consider  asymptotic  expansions  of  I  for  small  values  of  e. 


3.2.2  Asymptotic  Expansion  of  Minimum  Mean  Square  Error 
By  rewriting  the  numerator  in  Eq.  3-6  as 


[(l-e)bQ  <t>o  +  e  bi  4>i  ]2  -  [bQ[(l-e)  +  e  ft]  +  e(bi~bo)  ]2  , 


we  find  that 


X  -  b02  /  52[(l-e)$oU)  +  e  fc(0]d5 


+  2e  bo(b]-bo)  /  C2  $l(Ode 

C2  4>l(Q2 

+  e2(b!-b0)2  /  — - — — -  d5 

(l-e)<J>o(D  +  e  $!(£) 


(3-7) 


The  first  two  integrals  in  Eq.  3-7  can  be  computed  exactly;  the  third  cannot. 
Equation  3-5  and  3-7  give  us  the  following  asymptotic  expansion  of  the  minimum 
mean  square  error  V. 


o2  0Q2 


o4  [oi 2  -  on2  ] 


o4  [oi2  -  on2  ]2 


ALPHATECH,  INC. 


where 


“  C2  4>i(02 

J  =  /  -  d$ 

-co  (l-e)*0(O  +  e  4-1(0 


Since  the  third  term  of  Eq.  3-8  is  negative,  we  immediately  obtain  an  upper 


bound  of  V  which  we  denote  by  Vi+: 


o2  oq2 

Vj+  „  -  +  c 

o2  +  OQ2 


a4  [oi2  -  oo2  ] 

[O2  +  o02]2 


A  lower  bound  is  obtained  by  noting  that 


J  <  -  •  /  S2  4>l(Od5  -  e-l[o2  +  0l2]  . 

£  —  oo 


Thus,  the  lower  bound  is 


a2  0q2 


°M°12  ~  °02  1  /  °1  ”  °02 

+  e  -  I  1 - 


,2  +  0()2  [o2  4-  oq2]2 


+  a. 


Let  us  now  study  the  third  term  of  Eq.  3-8.  First  note  that  if 


oq2  <  o2  +  2oi^ 


then  the  integral 


°°  4>l(02 

/  C2  -  dc 

-co  4o(0 


diverges  and 


lira  J  = 


(3- 
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On  the  other  hand,  from  Eq.  3-11  it  is  clear  that 

lim  e  J  =  0  .  (3-14; 

e+0 

Equations  3-13  and  3-14  imply  that  the  higher  order  term  e2  j  in  Eq .  3-8 
decreases  more  slowly  than  e2  but  more  rapidly  than  e.  In  Appendix  A  we  show 
that  in  the  case  o  ■  »  oq,  the  correction  term  has  an  asymptotic  behavior 

like 

e2  J  ~  ]£n  e]1/2  e1+Y  (3-15] 


for  small  e,  where 

o2  +  Oq2 

Y  »  -  .  (3-16] 

2  2 
oiz  -  OQ* 


Appendix  A  also  gives  expressions  for  a  second-order  upper  bound  V2+  and  a 
second-order  lower  bound  V2~-  As  Appendix  A  makes  evident,  detailed  analytic 
approximation  of  V  can  be  difficult  for  even  very  simple  hybrid  state  estima¬ 
tion  problems. 


3.3  CRAMER- RAO  LOWER  BOUND 

Van  Trees’  [4]  extension  of  the  classical  Cramer-Rao  bound  provides 
another  method  of  approximating  the  minimum  mean  square  error  V.  This  lower 
bound,  which  we  denote  by  Vc“,  is  given  by 

vc-  -  (i  +  0-2  r1  (3-17; 


where  I  is  the  Fisher  information  given  by 


I  =■  E 


(3-18; 


-.yvvy V-V-V  V-  .  .  -V- 


V.' 


U  MU  I  BJ  J  HIUVIL  HJIW !  IV  l'*J  f  I AJ  !  IL'IVIW1  WI'W^'WWUW'l  VIW  l  V'V'f  *  1  ^  m*\  mj  |  ".'I' 
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and  p  is  the  conditional  density  of  the  measurement  y  given  the  state  x: 

p  *  p(y|x)  • 

In  this  simple  example,  p  is  the  Gaussian  sum  distribution 

p(y|x)  -  <l-e)^(5)  +  e  <!*(£)  (3-19) 


where 


£  =  y  -  x  , 

+i(£)  -  (2it  a±2)~  1/2  exp  {  -  1/2  Oi"2  £2}  . 


(3-20) 


It  is  easy  to  see  that  the  Fisher  information  is  given  exactly  by  the  integral 


I 


*  £2[(l-e)do“2  \(q(£)  +  e  oi~2  ^(£)  ]2 

/  -  d£ 

(l-e)'(O(C)  +  e  *i(£) 


(3-21) 


This  is  the  same  type  of  integral  that  occurred  in  the  exact  calculation  of  V. 
Although  we  cannot  express  this  integral  in  closed-form,  we  can  approximate  it 
in  the  same  way  we  approximated  V.  Thus,  we  find  that 

I  *  oo“2  -  e  oo“^(oi^  -  0Q^)  +  e2(°l~2  -  0q-2)2k  (3-22) 


where 


K 


r  &  wo2 

r  -  de  . 

-«  (l-e)'to(£)  +  e  i|q(£) 


Just  as  before  (Eqs.  3-9  and  3-11)  we  can  derive  an  upper  and  lower  bound 
of  K: 


0  <  K  <  e_1  Oi2  . 
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Thus,  we  obtain  the  corresponding  bounds  on  Vc“: 


Vr“  < 


o2  oq2 

32  +  an2 


1  -  e 


J*  [oi2  -  oo2  ] 

[a2  +  oq2  ]2 
a2  [ox2  -  oo2] 
oq2  [oq2  +  a2  ] 


(3-23) 


Vc"  > 


02  OQ2 
3 2  +  On2 


[°12  “  oo2]  oo2  <A 

j,  — -  i  ■—  •  ■ 

[ o2  +  oo2  ]2  oi2 
[oi2  -  oo2]  o2 
[a2  +  oq2  ]  ax2 


(3-24) 


From  the  inequality  Eq.  3-23  we  also  see  that 


o2  OQ2 


o2  +  oQ2 


o*[oi2  “  °02] 

+  e  -  +  Q(e) 

1«2  -  «„2 12 


(3-25) 


where  o(e)  is  an  error  term  for  which 


lim  e“l  o(e)  ■  0 
e-H) 


Thus,  we  see  that  the  first  order  in  e,  Vc"  agrees  with  Vx+,  the  upper  bound 
of  Eq.  3-10  and  also  the  first  order  expansion  of  V  (as  seen  from  Eq.  3-8). 
This  suggests  that  Vc~  might  be  an  accurate  estimate  of  V,  if  we  could  compute 
Vc~  exactly.  Note  that  theory  tells  us  that  Vx+  >  V  >  Vc“.  Thus,  Vx+  is,  in 
fact,  an  upper  bound  of  Vc“  and  it  is  tighter  than  the  upper  bound  in  Eq.  3-23 
as  one  can  easily  see. 


.  '■ 


22 


ALPHATECH,  INC. _ 

One  could  develop  an  approximation  of  the  integral  K  appearing  in  Eq. 

3-22  and  so  develop  second-order  upper  and  lower  bounds  of  I  and  thus  Vc-  in 
the  same  way  we  do  for  V  in  Appendix  A.  We  will  not  do  so,  but  let  us  note 
that  such  analytic  approximations  of  the  Fisher  information  in  Gauss-Markov 
hybrid  state  estimation  problems  are  difficult  to  obtain  even  for  simple 
problems.  In  Section  4  we  will  consider  the  Cramer-Rao  lower  bound  for  more 
general,  dynamic  hybrid  state  problems  and  find  that  there  are  good  reasons 
to  work  to  develop  better  methods  to  approximate  the  Fisher  information. 

3.4  RATE  DISTORTION  LOWER  BOUND 

The  last  method  we  will  consider  is  based  on  rate  distortion  theory  [16}, 
a  branch  of  Information  theory.  Rate  distortion  theory  allows  us  to  derive  an 
analytic,  closed-form  lower  bound  Vg“  of  the  minimum  mean  square  error  V.  The 
basic  approach  is  straightforward.  The  theory  [16]  tells  us  that  the  minimum 
mean  square  error  V  satisfies  the  inequality 

R(V)  <  I(y;x) 

where  R(V)  is  the  rate  distortion  function  of  the  state  x  and  I(y;x)  is  the 
mutual  information  between  x  and  the  measurement  y.  For  a  scalar  Gaussian 
state,  the  rate  distortion  is  simply 

R(v)  .  -L  . 

This  gives  the  lower  bound 

V  >  cj2  exp  {-  2  I(y;x)  } 

We  cannot  compute  the  mutual  information  I(y;x)  exactly,  but  we  can  find  a 
simple  upper  bound  I+  of  I(y;x)  and  thus  obtain  a  lower  bound  of  V,  namely 
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V  >  VR"  =  a2  exp{-  2  I+}  .  (3-26) 

To  find  an  upper  bound  of  mutual  information,  note  that* 

I(y;x)  =  h(y)  -  h(y|x)  (3-27) 

where  h(y)  is  the  differential  entropy  of  the  random  variable  y  and  h(y|x) 
is  the  conditional  differential  entropy  of  y  given  x.  The  measurement  y  is  a 
sum  of  the  state  x  and  the  independent  noise  v  which  has  a  probability  density 

p(n)  -  e  Pi(n)  +  (l-e)po(n) 

where 

Pi(n)  -  (2n  o±2)~  1^2  exp {-  1/2  n2  a±2 }  . 

Because  x  and  v  are  independent, 

h(y|x)  -  h(v)  . 


The  differential  entropy  h(v)  is  ,  by  definition, 

OP 

h(v)  -  /  p(y)  in  p(y)  dy  . 

•  CO 

The  function  f(p)  «  -p  In  p  is  concave  in  p,  and  application  of  Jenson’s 
inequality  yields  the  inequality 

h(v)  >  e  h(vi)  +  (l-e)h(vo) 


where  vi  has  density  pi .  Because  vj  is  Gaussian, 

h(v^)  “  -  in(2ne  o^2) 


*See  [15]  or  [16]  for  the  basic  definitions  and  results  of  information  theory. 
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Thus,  we  have  the  lower  bound 

h(y|x)  >  c  £n(2ire  o^2)  +  — —  (l-e)  £n(2ne  oo2)  .  (3-2 

Another  basic  result  of  information  theory  gives  us  an  upper  bound  on  the 
differential  entropy  of  a  random  variable  in  terms  of  its  variance.  In  par¬ 
ticular,  we  have 


h(y)  <  -  An(2ire  A) 

2 


where  A  is  the  variance  of  y,  given  by 


A  ■  ©2  +  e  0^2  +  (i-e)oQ2 


Substituting  this  expression  in  Eq.  3-29  and  using  it  together  with  Eq.  3-28 
in  Eq.  3-27,  we  obtain  the  upper  bound  I+ 

I(y;x)  <  1+  *  -  in(2xe[<j2  +  c  cri^  +  (l-e)oo2]) 

2 

-  -  e  An(2ne  ox^)  -  -  (1-E)in(2xe  oq2) 

2  2 

Using  this  expression  for  I*  in  Eq.  3-26  gives  us  the  rate  distortion  lower 


bound 


o2  ox 2  e  oq2(1-e) 
a2  +  e  01^  +  (l-c)oo^ 


Note  that  to  first  order  in  e,  the  bound  Vr“  is  given  by 


o2  0Q2 


02  0Q2 


VR-  =  -  +  e  - ;  <  2  in 

°2  +  °o2  [°o2  +  °2^  ( 


^ ^  [°o2  +  °2 1  “  [°12  "  °o2  ]  | +  °(  e> 


ALPHATECH,  INC. _ 

3.5  NUMERICAL  EXAMPLE 

Figure  3-1  shows  the  mean  square  error  predictions  made  by  the  different 
performance  analysis  methods  discussed  in  this  section.  The  statistical 
parameters  of  the  variables  in  Eq.  3-1  were  chosen  so  that 

x  »  0  , 

a2  «  100  , 

oi^  *  100  , 

oo^  “  1 

The  false  alarm  probability  e  **  P{q=l}  was  varied  from  0  to  .5. 

Figure  3-1  shows  performance  predictions  obtained  from  a  direct  first- 
order  expansion  (an  upper  bound  Vi+  in  Eq.  3-10  and  a  lower  bound  Vi“  in 
Eq.  3-12),  a  direct  second-order  expansion  (an  upper  bound  V2+  and  a  lower 
bound  V2~  in  Eq.  A-25  of  Appendix  A),  and  a  rate  distortion  bound  (a  lower 
bound  Vr~  in  Eq.  3-30).  Note  that  the  direct  first-order  expansion  Vi+  is 
equal  to  the  first-order  expansion  of  the  Cramer-Rao  lower  bound  Vc~.  Mean 
square  error  in  Fig.  3-1  is  normalized  by  dividing  by  the  initial  variance  o^. 

For  small  e  the  direct  second-order  expansion  should  give  tight  upper  and 
lower  bounds  of  mean  square  error.  As  can  be  seen  In  Fig.  3-1,  for  e  >  .3  the 
second-order  expansion  becomes  less  accurate  than  the  first-order  expansion. 

Note  that  the  rate  distortion  lower  bounds  for  this  example  is  less  accu¬ 
rate  than  the  other  performance  predictions.  This  appeared  to  be  true  in  all 
numerical  examples  we  considered. 


INITIAL 
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Figure  3-1.  Mean  Square  Performance  Bounds 
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3.6  REMARKS 

In  this  section  we  have  considered  several  performance  analysis  techniques 
applied  to  a  simple  static  Type  I  hybrid  state  estimation  problem.  Exact 
closed-form  evaluation  of  minimum  mean  square  error  is  impossible  to  achieve 
even  for  this  simple  problem.  Typically,  one  encounters  integrals  which  have 
as  integrands  the  ratio  of  two  Gaussian  sum  distributions.  These  integrals 
appear  in  the  exact  error,  the  Cramer-Rao  lower  bound,  and  in  the  second-order 
asymptotic  approximation  (Appendix  A).  Approximating  these  integrals  is  dif¬ 
ficult  even  for  the  simple  problem  of  this  section.  Generalization  to  dynamic 
problems  of  interest  seems  unlikely. 

On  the  other  hand,  the  rate  distortion  lower  bound  is  easy  to  compute, 
but  its  accuracy  seems  very  poor  compared  to  the  other  approaches  described 
in  this  section.  The  first-order  asymptotic  approximations  are  both  easy  to 
compute  and  reasonably  accurate.  Unfortunately,  this  type  expansion  does  not 
appear  to  generalize  to  dynamic  problems  of  interest. 

In  the  following  sections  we  will  consider  performance  analysis  techniques 
for  more  general  dynamic  hybrid  state  estimation  problems.  Section  4  discusses 
a  Cramer-Rao  approach  for  Type  I  problems,  Section  5  discusses  a  rate  distor¬ 
tion  approach  for  Type  II  problems,  and  Section  6  discusses  a  Monte-Carlo 
approach  that  potentially  applies  to  any  hybrid  state  problem. 
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SECTION  4 

CRAMER-RAO  LOWER  BOUND  FOR  TYPE  I  PROBLEMS 

4.1  INTRODUCTION 

Because  Type  I  problems  can  be  formulated  as  a  partially  observed  Markov 
process  <x(t),y(t)>  for  which  x(t)  is  continuous,  it  is  theoretically  possible 
to  apply  classical  performance  analysis  techniques  to  such  problems.  In  this 
section  we  will  examine  the  applicability  of  Cramer-Rao  methods  to  Type  I 
problems,  and  particularly  to  the  subclass  of  Type  I  problems  described  by 
Eq.  2-9. 

The  Cramer-Rao  method  as  presented  in  Van  Trees  [4]  provides  a  lower 
bound  of  the  estimation  error  covariance  (for  which  lower  is  interpreted  in 
the  sense  of  symmetric  matrix  inequalities)  in  Bayesian  estimation  problems 
in  which  both  state  and  measurement  are  modeled  as  random  variables.  Thus, 
the  Cramer-Rao  method  gives  us  an  estimate  of  mean  square  error.  Although  the 
bound  is  given  for  static  problems  in  [4],  it  has  been  extended  to  various  dy¬ 
namic  filtering  problems  in  [ 5 ] — [ 7 ]  (Bayesian  framework)  and  [8]  (non-Bayesian 
framework).  The  method  is  useful  because  it  gives  us  a  computable  performance 
measure  for  a  large  class  of  nonlinear  estimation  problems.  In  particular, 
the  computational  complexity  of  the  bound  for  the  filtering  problems  described 
in  ( 5 ] — ( 8 ]  is  comparable  to  that  required  to  compute  the  error  covariance  for 
a  linear  Gaussian  filter  and  is  proportional  to  the  number  of  time  periods 
under  investigation.  This  is  significantly  less  complexity  than  that  of  the 
optimal  filters  fot  most  nonlinear  problems. 
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We  will  show  here  that  the  Cramer-Rao  method  applies  similarly  to  Type  I 
hybrid  state  estimation  problems  described  by  Eq.2-9.  Actually,  we  show  more 
generally  that  this  method  applies  to  any  estimation  problem  for  which  the 
state  process  satisfies  Eq.  2-9a  and  the  measurement  process  y(t)  is  such 
that  all  y(t) ,  t-1,2,...,  are  conditionally  independent  and  identically  dis¬ 
tributed  given  the  state  process.  We  also  assume  that  the  conditional  dis¬ 
tribution  of  y(t)  depends  just  on  x(t)  and  that  this  distribution  p(y(t)|x(t)) 
is  a  continuously  differentiable  function  of  x(t)  which  vanishes  as  |x(t)| 
tends  to  infinity. 

Let  J(t)  denote  the  information  matrix  given  by 

J(t)  -  E{p(y(t)|x(t)}“2px(y(t)|x(t))Tpx(y(t)|x(t))}  (4-1) 

where  E{  }  denotes  taking  the  expectation  with  respect  to  the  joint  measure¬ 
ment  and  state  processes,  Px(y|x)  Is  the  partial  derivative  of  p(y|x)  with 
respect  to  x  (a  row  vector  if  x  is  a  column  vector),  and  pxT  denotes  the 
transpose  of  px  (a  column  vector  if  px  is  a  row  vector).  In  terms  of  J(t) 
the  Cramer-Rao  lower  bound  on  the  average  error  covariance  P(t)  is  given  by 
the  linear  filtering  update  equations,  namely 

P(t+1)  -  [(AP(t)AT  +  Q)_1  +  J(t)]-1  (4-2) 

where  Q  is  the  covariance  of  the  process  noise  w(t)  in  Eq.  2-9a.  We  will 
present  only  a  sketch  of  the  proof  here.  One  expresses  the  filtering  problem 
of  Eq.  2-9  as  a  large  dimensional  static  estimation  problem  of  estimating 
the  entire  state  trajectory  X  given  the  measurement  trajectory  Y.  Using  the 
static  result  of  [4},  note  that  the  Cramer-Rao  bound  for  this  problem  is 
equal  to  the  bound  (Indeed,  equal  to  the  exact  error  covariance)  for  a  linear 
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Gaussian  estimation  problem  for  which  Eq.2-9a  holds  and  Eq.  2-9b  is  replaced 
by  the  linear  equation 


y(t)  3  C(t)x(t)  +  v(t)  (4-3) 

in  which  C(t)  and  the  covariance  R(t)  of  v(t)  are  chosen  so  that 

J(t)  3  C(t)TR(t)-lc(t)  .  (4-4) 

This  is  always  possible  for  any  information  matrix  J(t)  (note  that  it  is  not 
necessary  to  compute  C(t)  or  R(t)).  It  follows  that  the  filter  estimation 
error  covariance  of  the  hybrid  state  problem  is  bounded  below  by  the  corre¬ 
sponding  error  covariance  for  the  linear  problem  given  by  Eqs.  2-9a  and  4-3, 
and  therefore,  the  bound  can  be  computed  according  to  Eq.  4-2. 

4.2  COMPUTATION  OF  INFORMATION  MATRIX  FOR  GAUSSIAN  MIXTURES 

From  above  we  see  that  applying  the  Cramer-Rao  approach  depends  on  com¬ 
puting  the  information  matrix  J(t)  in  Eq.  4-1  for  static  problems  of  the  form 
of  Eq.  2-9b,  namely 

y  -  C(q)x  +  D(q)v  ,  (4-5) 

where  x  and  v  are  Gaussian  distributed  and  q  takes  the  discrete  values 
1,2,...,N  with  probabilities  b^  for  each  value  q  3  k.  We  will  assume  without 
loss  of  generality  that  v  has  0  mean  and  an  identity  covariance  matrix.  Let 
Rk  denote  the  covariance  matrix  of  D(k)v,  that  is  R^  3  D(k)D(k)^.  In  this 
case,  the  conditional  probability  distribution  p(y|x)  needed  in  Eq.  4-1  is 

N 

p(y|x)  *  l  bkPk(y|x)  (4-6) 

k3l 
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where  pk  is  Che  Gaussian  density 


Pk(y|x)  =  (det  2*  Rk)-1^2 

exp(-l/2[y-C(k)]TRk-1ty-C(k)] ]  . 


The  corresponding  derivative  px(y|x)  is 


Px(y|x>  “  I  bk[y-C(k)]  Rk-1C(k)pk(y | x)  . 
k=l 


Define  J(x)  for  each  x  as 


(4-7) 


(4-8) 


r.", 

& 


i 
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» 
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J(x)  -  E|p(y|x)~2px(y|x)TpxCy|x) |x}  (4-9) 

where  the  expectation  E{  | x }  is  conditioned  on  x.  Thus,  the  information  matrix 
J  is  found  by  integrating  J(x)  with  respect  to  the  Gaussian  distributed  x. 

The  expectation  in  Eq.  4-9  involves  integrating  the  ratio  of  two  sums  of 
Gaussian  densities  given  by  substituting  Eq.  4-7  and  4-8.  In  general  it  is 
not  possible  to  do  this  integration  in  closed  form  (in  terms  of  elementary 
functions),  nor  is  J(x)  integrable  with  respect  to  the  Gaussian  density  of  x 
in  closed  form.  Thus,  one  must  resort  to  some  type  of  numerical  integration 
to  evaluate  J(x)  and  J.  Since  the  integrals  are  multidimensional  and  have 
infinite  domains  of  integration,  we  have  chosen  to  use  the  control  variate 
Monte-Carlo  method  of  approximating  integrals  [9], [10]. 

4.3  REMARKS  ON  CRAMER- RAO  APPROACH 

The  advantage  of  the  Cramer-Rao  approach  is  that  the  complexity  of  the 
computation  is  essentially  limited  to  the  complexity  of  a  static  problem 
depending  on  the  dimensions  of  the  state  vector  x(t),  the  measurement  vector 
y(t),  and  the  number  of  discrete  states  q(t).  Once  this  static  problem  is 
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solved  (namely  the  computation  of  J(t)  in  Eq.  4-1),  the  dynamic  problem  is 
easily  solved  by  the  familiar  recursion  in  Eq.  4-2. 

The  disadvantage  of  the  Cramer-Rao  approach  for  hybrid  state  performance 
analysis  is  that  the  information  matrix  J(t)  of  Eq.  4-1  may  be  difficult  to 
compute  even  for  small  dimensional  static  problems.  In  general  this  computa¬ 
tion  requires  integrating  a  ratio  of  two  sums  of  Gaussian  density  functions, 
and  thus,  the  integral  has  no  closed  form  expression  in  terms  of  elementary 
functions.  This  is  in  marked  contrast  with  the  situation  for  nonlinear  fil¬ 
tering  problems  in  Gaussian  white  noise  where  J(t)  can  be  computed  in  closed 
form  for  a  large  class  of  nonlinearities  [11], [12].  As  noted  above,  we  have 
employed  Monte-Carlo  methods  to  compute  these  integrals,  but  even  these  meth¬ 
ods  become  very  computationally  expensive  if  the  static  dimensions  of  the 
problem  are  large.  This  would  be  the  case  for  realistic  multiobject  tracking 
problems  with  many  targets  and  many  returns  per  measurement  period. 

The  Cramer-Rao  method  produces  a  lower  bound  to  the  minimum  mean  square 
estimation  error;  however,  this  lower  bound  may  be  much  lower  than  the  minimum 
mean  square  error  —  that  is,  the  bound  can  be  very  loose.  In  principle,  the 
bound  becomes  tight  as  the  measurement  covariance  tends  to  zero.  More  pre¬ 
cisely,  the  conditional  covariance  of  y  in  Eq.  4-5  given  x  must  be  small.  In 
practice,  one  must  determine  how  small  by  numerical  computation  and  comparison 
to  more  accurate  performance  measures. 

We  wish  to  mention  two  previous  applications  of  Cramer-Rao  methods  to 
hybrid  state  estimation  problems.  If  the  measurement  is  purely  discrete  or 
is  the  sum  of  a  continuous  and  discrete  component,  one  can  compute  Cramer-Rao 
type  bounds  fairly  easily  in  a  large  number  of  cases  [5], [13].  In  fact,  in 
this  case  it  is  possible  to  deal  with  a  discrete  process  q(t)  that  is  dependent 
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on  the  continuous  process  x(t).  These  computations  depend  crucially  on  the 
fact  that  the  discrete  state  enters  only  additively  in  the  measurement  equa¬ 
tion.  The  multiplicative  dependence  we  have  assumed  severely  complicates  the 
problem. 

The  work  reported  in  [14]  is  much  closer  to  the  problems  presented  here. 
That  work  studies  a  specific  hybrid  state  parameter  estimation  problem  using 
the  Cramer-Rao  method  for  unknown,  nonrandom  parameters  (see  [4]).  The 
approach  in  [14]  is  somewhat  different  as  the  state  x(t)  is  n*'?  random  as  it 
is  here,  and  the  approach  is  a  batch  procedure  rather  than  the  recursive  pro¬ 
cedure  presented  here.  Furthermore,  uniformly  distributed  continuous  input 
noise  variables  are  used  while  we  only  use  Gaussian  distributed  input  noise. 
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SECTION  5 

RATE  DISTORTION  LOWER  BOUNDS 
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5.1  RATE  DISTORTION  METHOD 

Information  theory  provides  another  approach  to  bounding  estimation  per¬ 
formance  based  on  Shannon's  theory  of  rate  distortion  [15], [16].  The  basic 
mathematical  idea  behind  that  approach  is  that  the  mean  square  error  D  of  any 
estimator  of  a  random  variable  X  given  measurement  Y  is  constrained  by  the 
inequality 

RX(D)  <  I(X;Y)  (5-1) 

where  Rx  is  the  rate  distortion  function  of  X  and  I(X;Y)  is  the  mutual  infor¬ 
mation  between  X  and  Y  [16]  .  Solving  Eq.  5-1  for  the  largest  D  satisfying 

the  inequality  gives  a  lower  bound  of  the  minimum  mean  square  error  of  esti¬ 
mating  X  given  the  measurement  Y.  In  a  filtering  problem  one  would  let  X  be 

the  current  state  component  xj_(t)  of  interest  and  let  Y  be  the  measurement 
sequence  <y(t) ,y(t-l) , . . . ,y(l)>  up  to  the  current  time.  This  method  has  been 
applied  to  nonlinear  filtering  problems  with  continuous  states  and  measure¬ 
ments  generated  by  Gaussian  noise  inputs  in  [  17 ] —  £ 20 ] .  We  will  apply  the 
method  to  Type  I  and  II  hybrid  state  estimation  problems  in  this  section. 


5.2  RATE  DISTORTION  BOUND  FOR  TYPE  II  PROBLEMS 


The  difficulty  of  applying  the  rate  distortion  method  depends  on  the 
difficulty  of  calculating  the  rate  distortion  function  Rx  and  the  mutual 
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information  I(X;Y)  in  Eq.  5-1.  In  general  (and  this  is  true  for  the  hybrid 
state  estimation  problem),  these  expressions  cannot  be  computed  in  closed 
form.  Therefore,  we  look  for  an  upper  bound  Iy  of  the  mutual  information  and 
a  lower  bound  of  the  rate  distortion  function.  The  Eq.  5-1  is  replaced  by 

RL(D)  <  Iu  .  (5-2) 

Solving  Eq.  5-2  for  the  largest  D  that  satisfies  the  inequality  gives  a  lower 
bound  of  the  minimum  mean  square  estimation  error  which  is  lower  than  the 
bound  given  by  Eq.  5-1. 

Let  us  first  find  a  lower  bound  for  the  rate  distortion  function. 
Shannon's  lower  bound  of  the  rate  distortion  function  Rx  is 

RX(D)  >  h(PX)  -  1/2  log(2ir  eD)  (5-3) 

where  h(px)  denotes  the  differential  entropy  of  the  probability  density  px 
of  X  [16].  To  avoid  notational  complexity  but  without  loss  of  generality 
(as  will  become  evident),  let  x(t)  be  a  scalar  process.  The  density  Px(t) 
is  a  Gaussian  mixture,  and  although  we  cannot  compute  h(px(t)),  we  can  find 
a  lower  bound  for  it,  namely 

h(px(c))  >  Ih(Px(t)|q)  P(q»t)  (5-4) 

q 

where  q  denotes  the  discrete  sequence  <q( t)  , . . .  ,q(l)> ,  p(q,t)  denotes  the 
probability  of  this  sequence  occurring,  and  px(t)|q  *s  the  Gaussian  density 
of  x(t)  conditioned  on  the  sequence  q.  This  gives  us  a  lower  bound  of  the 
rate  distortion  function.  For  reference  below  let  P(t|q)  denote  the  condi¬ 
tional  covariance  of  x(t)  given  the  sequence  q. 
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Next  we  find  an  upper  bound  for  the  mutual  information.  The  basic  rela¬ 
tion  is 


l(x( t) ;y(t) , . . . ,y(l) )  <  l[x(t),...,x(l),q(t),...,q(l);y(t),...,y(l)) 


(5-5) 


The  right  side  of  Eq.  5-5  is  given  by 


l(x(t) . x(l),q(t),...,q(l);y(t),...,y(l)) 

-  h(y(t) , . . . ,y(l) )  (5-6) 

“  h(y(t),...,y(l) | x(t),...,x(l),q(t ),..., q(l)) 


where  the  expression  h(Y|X)  denotes  the  conditional  differential  entropy  of  Y 
with  respect  to  X  [16].  Let  M(t|q)  denote  the  covariance  matrix  of  the  mea¬ 
surement  sequence  <y(t) , . . . ,y(l)>  conditioned  on  knowing  the  discrete  sequence 
q  »  <q( t) , . . . ,q(l)> .  Then  the  differential  entropy  is  bounded  above  by 

h(y(t),...,y(l))  <1/2  log [det (2xe  l  p(q,t)M(t|q)  )  ]  (5-7) 

q 

and  the  conditional  differential  entropy  is  given  exactly  by 


h(y(t),...,y(l)|x(t),...,x(l),q(t) . q(l)) 

=*  1/2  l  p(q,t)  log  [det  (2 neR(  1 1 q)  )  ] 

q 


(5-8) 


where  R(t|q)  denotes  the  covariance  matrix  of  the  noise  sequence  D(q(t))v(t) 
given  the  discrete  sequence  q(t).  Eqs.  5-7  and  5-8  give  us  an  upper  bound  of 
the  mutual  information. 

Using  the  lower  bound  of  the  rate  distortion  function  and  the  upper  bound 
of  t.ie  mutual  information,  we  obtain  the  following  lower  bound  of  the  mean 
square  estimation  error  of  estimating  x(t)  given  the  measurements  up  to 


time  t. 
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Dt  >  exp(Ft  +  Gt  -  Ht)  (5-9) 

where  Ft,  Gt ,  and  Ht  are  given  by 

Ft  -  l  P(q,t)log(P(t|q))  (5-10) 

q 

Gt  -  I  p(q,t)log(det[R(t|q)])  (5-11) 

q 

Ht  -  log(det[£  p(  t  ,q)M( 1 1 q) ]  )  (5-12) 

q 

Note  that  if  the  noise  random  variables  v(t)  are  independent,  standard 
Gaussian  random  vectors,  then  R(t|q)  is  the  block  diagonal  matrix 

R( 1 1 q)  -  diag {d (q(k) )  D(q(k))T}  (5-13) 

and 

t 

det  R(t|q)  **  n  det  |D(q(k)  )  D(q(k)  )? }  .  (5-14) 

k=l 

If  the  q(k)  are  independent  and  identically  distributed,  then 

Gt  -  t  •  E{log  det (D(q(l) )  D(q(l))T}  .  (5-15) 

5.3  REMARKS  ON  RATE  DISTORTION  METHOD 

Note  that  no  assumption  was  made  about  the  statistics  of  the  discrete 
state  process  q(t)  except  that  they  are  independent  of  the  continuous  states. 
Thus,  the  rate  distortion  bound  allows  us  to  treat  Type  II  problems  in  which 
q(t)  is  a  finite  state  Markov  chain.  In  theory,  one  could  treat  more  general 
problems  (Type  III)  using  the  rate  distortion  theory  approach  provided  that 
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one  can  compute  or  bound  the  corresponding  rate  distortion  function  for  the 
state  component  and  the  performance  criterion  of  interest.  For  example,  one 
can  do  this  if  one  is  interested  in  just  estimating  discrete  states  and  the 
average  probability  of  error  criterion  is  used.  See  [16]  for  the  general 
theory  of  rate  distortion  with  general  error  criteria. 

The  rate  distortion  method  has  a  different  kind  of  computational  complex¬ 
ity  from  the  Cramer-Rao  method  presented  above.  One  can  compute  this  bound 
exactly  (i.e.,  no  numerical  integrations  are  required),  but  the  complexity 
increases  exponentially  as  the  number  of  time  periods  t  considered.  The  com¬ 
putations  are  similar  in  nature  to  the  computations  occurring  in  the  optimal 
hybrid  state  estimation  algorithm  [1],[2]  —  that  is,  for  each  sequence  of 
discrete  states  one  computes  M(t|q)  and  P(t|q)  from  the  corresponding  linear 
filter  equations.  Of  course,  these  computations  are  exact  and  no  simulation 
or  Monte-Carlo  averaging  is  Involved.  Still,  if  there  are  N  discrete  values 
of  q(t),  then  this  computation  will  involve  Nc  filter  computations  to  consider 
a  time  period  of  t  steps.  Furthermore,  the  covariance  M(t|q)  has  dimension 
mt  where  m  is  the  dimension  of  y(t)  and  t  is  the  number  of  time  periods  con¬ 
sidered.  Computing  the  determinants  needed  in  Eqs.  5-11  and  5-12  requires 
on  the  order  of  (mt)^  operations  (see  [20], [21]  for  an  efficient  method  for 
doing  this  in  the  context  of  rate  distortion  problems).  Thus,  although  the 
rate  distortion  method  is  easier  than  the  Cramer-Rao  method  to  compute  for 
short  periods  of  time,  the  rate  distortion  computation  quickly  overtakes  the 
Cramer-Rao  in  complexity  as  the  number  of  time  periods  increases.  This  method 
has  a  similar  sensitivity  to  the  number  of  discrete  states  and  the  dimensions 


of  the  state  and  measurement  vectors. 
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The  rate  distortion  method  gives  a  lower  bound  of  minimum  mean  square 
error,  and  in  some  cases  that  lower  bound  may  be  very  loose.  In  principle, 
the  bound  becomes  tight  as  the  measurement  covariance  increases,  but  in  prac 
tice  (as  with  the  Cramer-Rao  method)  one  has  to  check  the  bound  with  more 
accurate  performance  measures  computed  for  some  simple  test  cases  in  order 
to  get  a  feel  for  the  rate  distortion  bound's  accuracy. 
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SECTION  6 

MONTE-CARLO  PERFORMANCE  ANALYSIS 
BASED  ON  THE  REPRESENTATION  THEOREM 


6.1  INTRODUCTION 

A  major  disadvantage  of  the  rate  distortion  method  described  in  Section  5 
is  that  it  can  be  inaccurate  unless  measurement  noise  is  very  high.  In  [28] 
we  faced  a  similar  problem  in  applying  rate  distortion  methods  to  nonlinear 
filtering  problems  such  as  occur  in  passive  tracking  problems.  In  this  sec¬ 
tion  we  will  extend  a  Monte-Carlo  method  discovered  in  [28]  to  include  hybrid 
state  estimation  problems.  Note  that  this  method  is  based  on  a  rate  distor¬ 
tion  approach  due  to  Caldos  [18], [19].  This  method  uses  the  representation 
theorem  of  nonlinear  filtering  theory  to  obtain  an  approximation  of  a  rate 
distortion  lower  bound  on  estimation  error.  Galdos  presented  both  a  pertur¬ 
bation  approximation  and  a  Monte-Carlo  approximation  for  his  method  in  [18] . 

In  [19]  he  clarified,  extended,  and  simplified  the  Monte-Carlo  approach.  In 
this  section  we  show  that  one  can  use  the  Monte-Carlo  approach  directly  to 
approximate  the  minimum  mean  square  estimation  error  without  invoking  any  rate 
distortion  theory.  As  we  will  see,  this  direct  method  not  only  avoids  the 
lower  bound  aspect  of  rate  distortion  theory,  but  is  also  is  simpler  to  imple¬ 
ment  and  more  generally  applicable  than  the  indirect  rate  distortion  approach 
of  Galdos.  In  subsection  6.2  we  will  generalize  Galdos'  method  to  general 
hybrid  state  systems.  In  subsection  6.3  we  will  present  our  direct  method, 
and  in  subsection  6. A  we  will  make  some  concluding  remarks. 


6.2  GALDOS  RATE  DISTORTION  METHOD 

The  problem  of  [19]  is  that  of  estimating  optimal  performance  in  discrete 
time  nonlinear  filtering  problems  described  by  equations  of  the  form 

x(t+l)  -  f(t,x(t))  +  g(t,x(t)  )w(t)  (6-1) 

y(t)  =  h(t,x(t) )  +  v(t)  (6-2) 

where  w(t)  and  v(t)  are  zero  mean  Gaussian  white  noise  sequences.  In  partic¬ 
ular,  if  d(x(t))  is  a  scalar  function  of  the  state  x(t),  it  is  desired  to 
approximate  the  minimum  mean  square  error  for  estimating  d(x(t))  given  mea¬ 
surements  y(s)  for  s  *  1 . t. 

The  method  of  [19]  is  not  limited  to  systems  of  the  form  of  Eqs.  6-1  and 
6-2.  As  we  will  show,  one  can  consider  any  problem  such  that  the  measurements 
y(s)  are  conditionally  independent  given  the  state  process  x(s),  l<s<t,  and 
the  conditional  density  that  y(s)  -  y  given  x(s)  -  x  is  p(y|x,s).  One  also 
needs  to  be  able  to  simulate  independent  samples  of  measurement  and  state 
processes  —  we  will  explain  more  precisely  what  is  needed  below. 

Rate  Distortion  Approach  [16], [19] 

One  can  consider  the  problem  of  estimating  the  minimum  expected  value 
of  p(d(x(t)) ,d(t) )  where  p  is  some  arbitrary  error  criterion.  The  minimum 
is  taken  over  all  estimates  d(t)  which  depend  only  on  the  measurements 
y(l),...,y(t).  The  rate  distortion  approach  is  based  on  the  inequality 

R(et)  <  l(d(x(t));Yc)  (6-3) 


where 
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R(D)  =  rate  distortion  function  with  respect  to  p; 

et  =  E{p(d(x(t)) ,d(t) ) }; 

d(t)  =  estimate  of  d(t)  based  on  Yt; 

Yt  3  measurement  sequence  y(l) , • . . ,y(t) ; 
p  =  error  criterion; 

I  =  mutual  information  between  d(x(t))  and  Yt . 

To  use  Eq.  6-3  you  must  be  able  to  compute  R  (or  rather  its  inverse  function 
R~l)  and  the  mutual  information  l(d(x( t) ) ; Yt  )•  Then 

et  >  R-1(l(d(x(t));Yt>)  . 

If  we  cannot  compute  the  rate  distortion  function  or  the  mutual  information, 
then  we  must  approximate  by  means  of  a  lower  bound  on  the  rate  distortion 
function, 

RL(D)  <  R(D) 

and  an  upper  bound  on  the  mutual  information. 

Iu(t)  >  I (d(x( t) ) ; Yt )  . 

This  gives  the  looser  inequality 

Et  >  RL-1(lu(t))  (6-4) 

for  the  estimation  error. 

Galdos  uses  the  usual  lower  bound  R^  on  the  mean  square  distortion  func¬ 
tion  R  for  the  error  criterion  p(d,d)  3  (d-d)2  and  the  source  d(x(t))  given  by 

RL(D)  =*  h(pd(x(t))) - —  log(2neD)  (6-5) 
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where  Pd(x(t))  denotes  the  probability  density  of  d(x(t))  and  h(p)  denotes  the 
differential  entropy  of  a  probability  density  p.  Note  that  Rl  (and  R  also) 
depends  on  t  implicitly  through  this  differential  entropy  of  d(x(t)).  Note 
that  if  d  is  a  linear  functional  of  x(t)  and  if  x(t)  is  Gaussian  distributed 
(as  it  is  for  the  subclass  of  Type  I  hybrid  state  systems  described  in  Eq.  2-9, 
for  example),  then  the  differential  entropy  is  easily  computed  as  a  function 
of  the  variance  V  of  d(x(t)),  namely 


h(Pd(x(t)))  -  ~  log(2weV)  . 


(6-6) 


The  mutual  information  is  not  so  easy  to  approximate,  and  this  is  where 
Galdos  resorts  to  Monte-Carlo  approximation.  By  definition,  the  mutual  infor¬ 


mation  is 


l(d(x(t));Yt]  *  E<  log 


Py |d  lYt Id (x(t) ) } 
Py  {Yt  1 


(6-7) 


In  the  right-hand  side  expression  of  Eq.  3-7,  we  use  Py|d  to  denote  the  con¬ 
ditional  density  of  the  measurement  sequence  Yt  given  d(x(t)),  and  we  use  py 
to  denote  the  unconditional  density  of  Yf 

Representation  Theorem 

The  representation  theorem  (for  discrete-time  estimation,  this  is  just 
Bayes'  rule)  is  used  to  represent  the  ratio 


Py|d{*tid(x(t))}  Ex|d{exp(ct)|d(x(t) )} 
Py  (Yt }  Ex{expUt)} 


(6-8) 


where  is  the  random  variable  defined  by 
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t 

^t(YfXt)  -  I  log[p(y(s)|x(s),s)]  .  (6-9) 

s=*  1 

Thus,  Ct  is  a  function  of  the  measurement  sequence  Y^  and  the  state  sequence 
Xt  *  (xi,...,xt).  The  expectations  Ex|d  and  Ex  integrate  over  the  state 
sequence  only,  leaving  the  measurement  sequence  alone.  To  be  precise,  let  px 
denote  the  probability  density  of  the  state  sequence  Xt .  Then  Ex  is  inter- 
pre.ed  to  mean 

Ex{expUt)}  -  /  exp(ct(Yt,Xt))  px  {Xt }  dXt  .  (6-10) 

Thus,  this  expectation  produces  a  function  of  Yf  If  Px|d  denotes  the  condi¬ 
tional  density  of  the  state  sequence  given  the  variable  d(x(t)),  then  Ex|d  is 
interpreted  to  mean 

Ex|d{exp(i;t)}  -  /  exp(ct0ft,Xt))  Px|d  {Xt  N  (x(  t)  ) }  dXt  .  (6-11) 

This  expectation  produces  a  function  of  Yj;  and  d(x(t)). 

Indirect  Monte-Carlo  Approximation 

The  approach  of  [19]  is  to  approximate  the  integrals  in  Eqs.  6-10  and 
6-11  by  Monte-Carlo  methods,  and  then  use  Monte-Carlo  methods  again  on  the 
approximate  ratio  corresponding  to  Eq.  6-8  to  obtain  the  final  approximation 
for  the  mutual  information  in  Eq.  6-7.  We  will  now  outline  a  simple  Monte- 
Carlo  algorithm  similar  to  [19]  for  approximating  the  mutual  Information. 

Note  that  the  integers  n  and  m  control  the  number  of  simiulations  used  in  the 
Monte-Carlo  approximation. 

Step  1.  Initialize  Ic  0  and  i  ♦  1 . 

Step  2.  Simulate  d(x(t)j  and  the  measurement  sequence  y( 1 ) , . . . ,y( t) 
from  the  joint  state  and  measurement  model. 
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Step 

3. 

Initialize  Nt  ♦  0,  Dt 

Step 

4. 

Simulate 

the  state  sei 

d(x' ( t) ] 

=  d{x(t) ). 

j 

1 

Step 

5. 

Nt  +  — 

Nr  +  -  exp 

j+1 

j+1 

Step 

6. 

Simulate 

a  new  state 

conditioning . 

j 

1 

Step 

7. 

Dt  ♦  - 

Dt  +  -  exp 

j+1 

j+1 

Step 

8. 

If  j  -  «, 

then  go  to 

Step  4. 

i 

1 

Step 

9. 

It  *  - 

It  +  -  log 

i+1 

i+1 

Step 

10. 

If  i  =*  n, 

then  stop. 

When  the  algorithm  stops,  the  value  of  It  will  be  an  approximation  of  the 
mutual  information  in  Eq.  6-7.  Using  this  with  Eq.  6-5  in  the  inequality  of 
Eq.  6-4  gives  the  approximate  lower  bound 


et  >  exp{2  [h(pd(x(t))  )  "  It]}  (6-12) 

i  ite 


for  the  mean  square  error  ec  of  estimating  d(x(t)). 

Remarks  on  Indirect  Monte-Carlo  Approach 

A  few  remarks  are  in  order  at  this  point.  The  algorithm  consists  of  an 
outer  loop  that  iterates  n  times  (Steps  2  through  10)  and  an  inner  loop  that 
iterates  m  times  (Steps  4  through  8)  for  each  iteration  of  the  outer  loop. 
Thus  the  algorithm  requires  of  the  order  of  nm  simulations  (it  requires  n 
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unconditional  simulations  of  the  state  and  measurement  sequences,  nm  uncondi¬ 
tional  simulations  of  the  state  sequence,  and  nm  conditional  simulations  of 
the  state  sequence).  As  n,m  *  °»,  the  approximation  It  converges  to  the  mutual 
information.  However,  the  rate  of  convergence  is  not  clear.  For  some  numer¬ 
ical  experiments,  n  and  m  must  approach  1000  for  reasonable  accuracy.  Thus, 
the  simple  Monte-Carlo  approach  outlined  above  can  be  very  time-consuming, 
although  the  algorithm  is  very  simple  to  implement. 

A  more  difficult  problem  is  simulating  the  conditional  state  sequence  in 
Step  4.  Theoretically,  this  involves  simulating  a  time  reversed  state  x(t) 
process  [19], [29]  from  a  terminal  condition  specified  on  d(x(t)).  If  the 
state  process  is  Gaussian  and  d  is  a  linear  functional  of  the  state  vector, 
then  the  reversed  process  can  be  simulated  by  a  linear  Gaussian  model  as  out¬ 
lined  in  [19]  and  derived  in  [29.  If  the  state  is  non-Gaussian  or  nonlinear, 
or  d  is  nonlinear,  then  the  backward  simulation  problem  is  considerably  more 
difficult.  Galdos  suggests  an  approximate  approach  in  [19]. 

Computing  the  differential  entropy  h(pd(x(t)))  easily  accomplished  by 
using  Eq.  6-6  if  the  state  process  is  Gaussian  and  d  is  linear.  Note  also 
that  in  the  linear  case  the  approximation  Rl  in  Eq.  6-5  is  exact:  Rj,  •  R. 

If  one  or  the  other  is  nonlinear,  then  the  differential  entropy  may  also  have 
to  be  approximated. 

Finally,  the  method  is  fundamentally  limited  by  the  fact  that,  at  best, 
it  is  a  lower  bound  on  performance.  There  is  no  guarantee  that  the  rate  dis¬ 
tortion  lower  bound  will  be  a  good  approximation  of  the  estimation  error. 

Despite  these  limitations,  the  method  has  some  appeal  because  it  poten¬ 
tially  applies  to  a  very  general  class  of  filtering  problems  and  requires 
little  special  analysis  once  the  basic  algorithm  is  set  up.  In  the  next 
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section  we  will  show  how  to  use  the  representation  theorem  in  a  direct  Monte- 
Carlo  approach  that  avoids  many  of  the  problems  inherent  in  the  indirect  rate 
distortion  approach. 


6.3  DIRECT  MONTE-CARLO  PERFORMANCE  ANALYSIS 

The  direct  method  we  propose  is  based  on  the  observation  that  the  repre¬ 
sentation  theorem  used  in  [19]  essentially  gives  the  conditional  probability 
density  i.t  of  x(t),  and  this  density  allows  one  to  compute  the  exact  minimum 
mean  square  error  et  in  estimating  d(x(t))  via  the  formula 

et  =  E  j  /  d(C)2*t(Od5  -  (/  d(D*t(Od02  j  •  (6-13) 


To  see  how  this  is  done,  consider  the  representation  theorem  of  Eq.  6-8  again. 
Let  Px(t)  denote  the  unconditional  density  of  x(t)  and  let  Ex|x(t)  denote  the 
conditional  integral  defined  in  Eq.  6-11.  Then  the  conditional  density  itt  is 
given  by 

Ex|x(t){expUt)Ul 

"tU) '  ~i«>ut)i  *(c>tu  •  <6-u> 


Thus,  the  conditional  integrals  in  Eq.  6-13  can  be  expressed  as 


/  d(0*tU)dC  =  J  d(C) 


Ex|x(t)  {exp( Ct) I  5 } 
Ex{exPUt)  } 


Px(t)U)d5  , 


(6-15) 


which  is  same  as 


/  d(C)irt(OdC 


Ex{d(x(t)  )exp(ct)  } 
Ex  (expUt)  } 


(6-16) 


i 


Similarly,  we  have 
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J  dU)2"t(Ode  = 


Ex{d(x(t)  )2exp(ct)  } 

EX{exp(ct)  } 


We  will  now  use  these  formulas  to  derive  a  direct  Monte-Carlo  approach  to 
approximating  minimum  mean  square  error. 


Direct  Monte-Carlo  Algorithm 

We  describe  a  simple  Monte-Carlo  algorithm  analogous  to  the  one  given 
for  the  rate  distortion  approach.  As  above,  the  integers  n  and  m  determine 
the  number  of  simulations  used  by  the  Monte-Carlo  approximation. 

Step  1.  Initialize  ec  0,  i  ♦  1. 

Step  2.  Simulate  a  measurement  sequence  y(s),  l<s<t. 

Step  3.  Initialize  Vt  ♦  0,  Mt  ♦  0,  Pt  ♦  0,  j  «-  1. 

Step  4.  Simulate  the  state  sequence  x(s),  l<s<t. 
t 

Step  5.  ♦  I  log[p(y(s)|x(s),s)]  . 

s=-l 

Step  6.  Vt  ♦  — Vt  +  -  d(x(t) )2  exp(ct)  . 

j+1  j+1 

j  1 

Step  7.  Mt  ♦  -  Mt  +  d(x(t) j  exp(^t)  • 

j+1  j+1 

j  1 

Step  8.  Pt  ♦  Pt  +  exp(ct)  . 

j+1  j+1 

Step  9.  If  j  ■  m,  then  go  to  next  step,  else  j  ♦  j+1  and  return  to 
Step  4. 

i  1  pt  f  «t  ~12“ 

Step  10.  et  -  -  er  +  -  —  -  — 

i+1  i+1  [_Pt  L  pt  J 

Step  11.  If  i  *  n,  then  stop,  else  1  «-  i+1  and  return  to  Step  2. 


The  final  value  of  et  will  be  an  approximation  to  the  minimum  mean  square 
filtering  error  for  estimating  dfx(t)). 
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Remarks  on  Direct  Monte-Carlo  Approach 

The  direct  Monte-Carlo  approach  outlined  above  also  requires  of  the  order 
of  nm  simulations  (n  simulations  of  the  state  and  measurement  sequences,  and 
nm  simulations  of  the  state  sequence  alone).  Furthermore,  the  direct  method 
can  converge  slowly  just  as  the  indirect  method,  and  n  and  m  may  have  to  be 
very  large  to  obtain  sufficient  accuracy. 

However,  in  other  respects,  the  direct  method  is  an  improvement  over  the 
indirect  method.  The  direct  does  not  require  conditional  simulation  of  the 
state  sequence;  only  forward  time  simulations  are  needed,  nor  does  the  method 
require  any  entropy  calculations.  Consequently,  the  direct  method  is  easy  to 
apply  to  problems  with  non-Gaussian  noise  and  nonlinear  state  dynamics.  In 
fact,  the  type  of  state  model  is  irrelevant  to  the  application  of  the  algo¬ 
rithm,  provided  one  can  simulate  the  state  sequence  and  provided  that  the  mea¬ 
surement  model  is  given  the  probability  density  function  p(y|x,s).  Thus,  it 
is  possible  to  apply  the  direct  method  to  very  general  non-Gaussian,  non- 
Markovian  state  processes. 

Finally,  the  direct  method  approximates  the  actual  minimum  mean  square 
error  and  not  a  lower  bound  of  this  error.  Thus,  as  the  Monte-Carlo  parameters 
n  and  m  increase,  one  obtains  better  approximations  of  the  optimal  performance. 
For  the  indirect  method,  the  performance  approximation  may  be  poor  no  matter 
how  large  n  and  m  are . 

6.4  CONCLUSIONS 

In  [18], [19]  Galdos  had  the  innovative  idea  of  using  the  representation 
theorem  of  nonlinear  filtering  theory  for  computational  rather  than  theoret¬ 
ical  purposes  to  obtain  a  rate  distortion  lower  bound  of  optimal  filtering 
error  approximations  for  a  large  class  of  filtering  problems.  In  this  section 
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we  have  shown  how  to  apply  the  representation  theorem  directly  without  having 
to  use  rate  distortion  theory,  thus  obtaining  a  simpler,  more  accurate  approx¬ 
imation  of  optimal  filtering  error  which  applies  more  easily  to  a  larger  class 
of  hybrid  state  filtering  problems. 

The  technique  outlined  in  the  previous  subsection  has  the  advantage  that 
it  is  very  to  use  for  almost  any  filtering  problem.  Its  major  disadvan¬ 

tage  is  due  to  the  slow  convergence  of  the  crude  Monte-Carlo  methods  employed. 
This  suggests  a  problem  for  further  work.  The  paper  [30]  of  LeGrand  is  of 
interest  in  studying  this  problem  as  he  has  investigated  the  problem  of  accel¬ 
erating  Monte-Carlo  convergence  for  conditional  expectation  type  calculations 
using  importance  sampling.  Note  also  that  the  algorithm  described  in  subsec¬ 
tion  6.3  can  be  partly  decomposed  into  parallel  computations,  so  the  algorithm 
is  a  good  candidate  for  array  processing. 

Any  simulation  technique,  such  as  the  one  described  in  this  section,  has 
the  fundamental  limitation  that  it  only  indicates  performance  for  a  single  set 
of  parameters  defining  the  statistics  of  the  underlying  model.  Thus,  such 
methods  cannot  easily  provide  insight  into  how  performance  is  affected  by 
parametric  variations.  Unfortunately,  accurate  analytic  methods  which  might 
provide  such  insight  are  lacking  for  most  non-Gaussian,  nonlinear  filtering 
problems.  The  method  discussed  here  offers  a  viable  approach  to  accurate 
optimal  performance  analysis  which  may  at  least  help  determine  when  practical 
nonlinear  filters  are  near-optimal.  If  the  computational  problems  described 
above  can  be  solved  by  algorithmic  or  computer  processing  improvements,  it  is 
also  possible  that  the  technique  might  give  useful  insights  into  hybrid  state 
estimation  comparable  to  what  is  possible  with  analytic  techniques  available 
for  simpler  problems. 
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SECTION  7 

A  RANDOM  POINT  PROCESS  APPROACH  TO  MULTIOBJECT  TRACKING 

7.1  INTRODUCTION 

This  section  presents  a  new  theoretical  approach  to  multiobject  tracking 
problems  which  might  provide  a  framework  for  developing  effective  analytical 
tools  for  multiobject  tracking  problems.  The  approach  is  based  on  a  random 
spatial  point  process  model  of  measurements  and  a  Laplace  transform  method  for 
calculating  the  likelihood  formulas  that  are  essential  to  deriving  tracking 
algorithms  and  to  analyzing  tracking  performance.  The  model  Includes  the  one 
described  above,  but  it  is  much  more  compact  to  express  and  simpler  to  under¬ 
stand  (once  easy,  but  unfamiliar  mathematics  is  grasped).  Moreover,  we  show 
how  to  use  simple  Laplace  transform  arguments  to  analyze  this  model  and  avoid 
some  of  the  notational  headache  inherent  in  previous  approaches. 

Theoretically,  our  approach  is  close  to  and  was  inspired  by  the  random 
space-time  process  work  of  [22],  [23].  However,  the  space-time  processes  and 
the  problem  described  in  this  section  are  different  from  [22], [23].  This 
section  considers  only  spatial  point  processes  and  time  is  discrete. 

In  the  next  subsection  we  present  the  mathematical  model  and  background 
concerning  random  point  processes  needed  to  understand  and  analyze  the  model. 
In  subsection  7.3  we  state  and  derive  the  likelihood  ratio  formulas  for  such 
random  point  process  measurements  using  a  simple  Laplace  transform  approach. 
Subsections  7.4  and  7.5  illustrate  the  use  of  the  resulting  likelihood  ratio 
formula  by  deriving  the  (optimal  Bayesian)  hypothesis  tree  type  algorithm 
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(subsection  7.4)  and  showing  how  to  compute  Cramer-Rao  lower  bounds  of  mean 
square  estimation  error  covariances  (subsection  7.5).  Subsection  7.6  con¬ 
cludes  the  section. 

7.2  RANDOM  POINT  PROCESS  MODELS 

Conceptually,  the  approach  of  this  section  is  to  consider  the  observa¬ 
tions  occurring  in  one  time  period  or  scan  of  data  as  an  image  of  points 
rather  than  as  a  list  of  measurements.  The  approach  is  intuitively  appealing 
because  it  corresponds  to  one's  natural  idea  of  radar  and  sonar  displays  — 
devices  that  provide  two-dimensional  images  of  each  scan  of  data.  We  will 
show  in  this  section  and  the  next  that  the  approach  is  also  technically 
appealing  because  it  allows  more  compact  and  clearer  formulation  of  models 
and  derivation  of  estimators.  The  approach  uses  the  basic  ideas  of  abstract 
measure  theory  but  only  the  simplest  results.  The  reader  can  find  the  neces¬ 
sary  theory  in  [24], [25]. 

Following  [22], [23],  we  formulate  the  image  mathematically  in  terms  of 
a  random  point  process.  The  precise  way  to  do  this  is  to  use  random  measures 
[22] , [23] , [26] .  It  is  mathematically  convenient  to  represent  an  image  of 
points  in  terms  of  a  measure  taking  nonnegative  integer  values.  If  p  is  such 
a  measure  defined  with  respect  to  the  measure  space  Y,  then  p(B)  is  the  num¬ 
ber  of  points  falling  inside  a  measurable  subset  B  of  Y.  We  will  call  such 
integer-valued  measures  counting  measures.  An  important  example  of  such  a 
measure  is  the  Dirac  delta  measure  denoted  6y  for  a  given  point  y  of  Y  and 
defined  so  that  <5y(B)  »  0  if  y  is  not  in  B  and  6y(B)  »  1  if  y  is  in  B.  If 
y(k),  l<k<n,  is  a  sequence  of  points  in  Y,  then  the  sum 

n 

u  ■  I  <5y(k)  (7-1) 

k=l 
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is  also  a  counting  measure.  We  will  use  such  sums  to  represent  all  the 
counting  measures  we  need  to  deal  with.  Note  that  the  representation  of  y 
by  the  sequence  y(k)  is  unique  except  for  arbitrary  ordering  of  the  sequence 
of  points  y(k). 

To  formulate  a  statistical  measurement  model  we  must  define  random 
counting  measures.  Random  counting  measures  can  be  defined  precisely  as 
random  variables  taking  values  in  M(Y),  the  collection  of  counting  measures 
defined  with  respect  to  the  measure  space  Y  [26].  In  particular,  if  y  is  a 
random  variable  with  respect  to  Y  then  Sy  defines  a  random  counting  measure 
with  respect  to  M(Y).  One  can  generate  other  random  counting  measures  by 
summing  several  such  Dirac  delta  measures  as  in  Eq.  7-1. 

The  most  familiar  and  important  example  of  a  random  counting  measures 
is  the  Poisson  measure.  This  is  called  a  Poisson  process  in  [26]  because  it 
generalizes  the  conventional  Poisson  process  defined  on  the  real  line.  The 
Poisson  measure  can  be  defined  generally  as  follows.  Let  X  be  a  finite  (non¬ 
negative,  real-valued)  measure  with  respect  to  Y.  Let  N  denote  a  Poisson 
random  variable  with  mean  value  equal  to  X(Y),  and  let  y(k),  l<k,  denote  an 
infinite  sequence  of  independent  random  variables  with  re  .  :ct  to  Y  which  are 
independent  of  N  and  which  all  have  the  probability  distribution  given  by 
X(Y)”^X.  The  summation 

N 

v  3  l  5y(k)  (7-2) 

k»l 

defines  a  Poisson  measure  with  Intensity  measure  X.  The  conventional  Poisson 
process  N(t)  is  related  by  n(t)  ■  v([0,t])  to  the  Poisson  measure  v  whose 
intensity  measure  is  the  Lebesgue  measure  on  the  real  line. 
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A  Poisson  measure  has  the  property  that  for  each  measurable  subset  B  of 
Y ,  the  integer-valued  random  variable  v(B)  is  Poisson  distributed  with  mean 
value  X(B).  In  particular,  one  has  that 

E { v( B )  }  =*  X(B)  .  (7-3) 

Note  that  a  Poisson  measure  also  has  the  property  that  the  random  variables 
v(B^)  and  v(B2)  are  independent  if  the  sets  B^  and  B2  are  disjoint.  These 
properties  make  the  Poisson  measures  a  reasonable  model  for  clutter  in  multi¬ 
object  tracking  problems.  The  Poisson  distributed  random  variable  v(B)  models 
the  number  of  false  alarms  appearing  in  the  subset  B  of  the  observation  space 
in  one  scan  of  data.  The  intensity  X(B)  is  thus  the  expected  number  of  false 
alarms  appearing  in  B  in  one  scan. 

We  can  now  specify  the  random  point  process  model  we  use  for  multiobject 
tracking.  Let  x(t)  denote  the  joint  vector  of  individual  target  states  x^(t), 
l<i<n,.  For  simplicity  we  will  assume  that  there  is  a  known  finite  number  n 
of  targets  and  that  the  joint  process  x(t)  is  a  Markov  chain  taking  values  in 
the  measure  space  X.  In  subsection  7.4  we  will  further  assume  that  the  chains 
x^(t)  are  independent,  but  no  such  assumption  is  needed  in  subsection  7.3. 

The  exact  nature  of  these  chains  does  not  affect  the  results  obtained  here; 
for  example,  one  can  suppose  that  they  are  independent,  Gaussian  processes 
generated  by  finite  dimensional  linear,  Gaussian  systems. 

The  measurement  process  is  a  sequence  pt  of  random  counting  measures 
which  are  conditionally  independent  given  the  state  process  x(t)  and  whose 
conditional  distribution  at  time  t  depends  only  on  x(t).  Define  pt  as  the  sum 
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where  vt  is  a  Poisson  measure  with  intensity  X,  and  Tt  is  a  random  counting 
measure  of  the  form 

n 

*t  “  l  zi<0  (t)  *  (7-5) 

i-1  1 

In  Eq.  7-5  the  random  variable  z^(t)  is  1  or  0,  depending  on  whether  the  ic^ 
target  is  observed  (i.e.,  has  a  return)  in  scan  t  or  not.  The  random  variable 
yji(t)  takes  values  in  Y  and  is  the  observation  one  receives  in  scan  t  from  the 
ich  target  provided  that  it  is  detected.  We  assume  that  the  Poisson  measures 
vt  are  Independent  of  each  other  and  of  the  processes  x(t),  z^(t),  and  y^(t) 
for  all  i  and  t.  The  random  variables  z^(t)  and  y^(t)  are  conditionally  inde¬ 
pendent  for  all  i  and  t  given  the  state  process  x(t).  Furthermore,  the  con¬ 
ditional  distribution  of  yi(t)  is  given  by  the  measure  B  -*■  Hi(B|x(t))  with 
respect  to  the  measure  space  Y,  and  the  conditional  probability  that  z^(t)  ■  1 
is  given  by  G^(x(t)].  We  will  make  further  assumptions  concerning  these  dis¬ 
tributions  in  subsections  7.3  and  7.4. 

The  model  in  Eq.  7-4  succinctly  represents  a  scan  of  data  that  consists 
of  false  alarms  (modeled  by  vt)  and  target  detections  (modeled  by  T^) •  The 
false  alarms  are  modeled  appropriately  by  the  Poisson  measure  vt  which  assumes 
as  a  reasonable  approximation  that  false  alarms  are  spatially  independent  of 
each  other.  Note  that  the  intensity  measure  X  allows  us  to  vary  the  expected 
rate  of  false  alarms  over  the  observation  space  Y.  For  example,  the  observa¬ 
tion  space  Y  might  simply  be  the  plane  and  X  might  be  a  constant  multiple  of 
the  uniform  distribution  over  a  disk  representing  the  area  of  coverage  of  the 
sensor.  This  implies  a  constant  rate  and  density  of  false  alarms  over  the 
sensor's  area  of  coverage.  The  target  detections  are  modeled  by  the  random 
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counting  measure  Tt  given  in  Eq.  7-5.  This  model  assumes  that  the  group  of 
targets  produces,  at  most,  n  measurements.  In  subsection  7. A  we  will  assume 
additionally  that 

Hi(B|x)  =  H(B|Xi)  (7-6) 

and 

Gt(x)  =  G(xi)  (7-7) 

so  that  each  target  produces,  at  most,  one  observation  and  each  observation 
comes  from,  at  most,  one  target.  The  probability  distribution  H(B|x^)  corre¬ 
sponds  to  the  observation  of  a  detected  target;  for  example,  it  might  be  a 
multivariate  Gaussian  distribution  whose  covariance  is  constant  and  whose  mean 
is  a  linear  function  of  x^(t).  The  quantity  G(x^)  is  the  detection  proba¬ 
bility  for  a  target  in  state  x^;  for  example,  this  might  be  a  constant  inde¬ 
pendent  of  X£. 

In  the  next  section  we  will  derive  a  likelihood  function  for  the  measure¬ 
ment  pt  given  the  state  x(t).  Before  doing  so,  we  need  to  present  a  few  more 
basic  mathematical  results  about  random  counting  measures.  Following  [26],  it 
is  convenient  to  denote  by  p(<j>)  the  ordinary  Lebesgue-Stielt jes  integral  of  a 
real-valued,  integrable  function  $  with  respect  to  a  measure  p.  That  is,  we 
define 

b(ij>)  a  dp  .  (7-8) 

If  p  is  a  counting  measure  given  by  Eq.  7-1,  then  this  integral  is  just  the 
finite  sum 

n 

U(<)>)  “  I  $(y(k))  •  (7-9) 

k=l 
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If  the  function  $  is  nonnegative,  the  integral  y(  $)  is  always  well-defined, 
although  it  may  be  infinite.  If  $  is  a  bounded  real  function  and  y  is  a 
finite  measure,  then  41  is  also  integrable.  If  y  is  a  random  counting  measure, 
a  stochastic  integral  can  be  defined  just  as  in  Eq.  7-9.  In  our  examples,  all 
random  counting  measures  are  finite  so  that  y(  $)  is  a  well-defined  integral 
for  bounded  real  functions  <p. 

Stochastic  integrals  such  as  y( 4)  are  important  because  we  use  them  to 
define  the  Laplace  transform  of  the  random  measure  y  by 

Ly(4>)  -  E {exp (- y(  <p)  ) }  (7-10) 

for  all  nonnegative  real  function  $  [26]  •  The  Laplace  transform  Ly  uniquely 
determines  the  distribution  of  the  random  measure  y  and  is  often  easier  to 
work  than  the  probability  distribution  of  a  random  measure.  An  important 
special  case  is  the  Laplace  transform  of  the  Poisson  measure  v  given  by 

LVU)  -  exp(-X(l-exp(-4>))  )  (7-11) 

where  X  is  the  intensity  measure  of  v  and  X(l-exp(-$) )  is  the  integral  of 
1  -  exp(-iji)  with  respect  to  X.  Note  that  if  X  is  a  finite  measure,  then  v  is 
almost  surely  a  finite  counting  measure  and  Eq.  7-11  is  true  for  bounded  <f>  as 
well  as  nonnegative  $.  As  a  final  example,  suppose  that  t  Is  the  Dirac  delta 
measure  5y  and  y  is  a  random  variable  with  probability  measure  H.  Then  the 
Laplace  transform  of  t  is 

Lt(4>)  -  H(exp(-$))  .  (7-12) 

In  the  next  subsection  we  will  use  the  results  (Eqs.  7-11  and  7-12)  to  derive 
likelihood  function  for  the  multiobject  tracking  measurement  model  as  expressed 
In  Eqs.  7-4  and  7-5. 
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7.3  LIKELIHOOD  RATIO  FUNCTIONS 

The  key  to  developing  an  estimation  algorithm  for  the  measurement  model 
given  in  Eqs.  7-4  and  7-5,  is  to  express  the  likelihood  p(ptlx(t))  of  observing 
given  the  true  state  is  x(t).  Since  yt  is  a  rather  complicated  mathematical 
object  (namely  a  counting  measure),  it  is  not  obvious  what  the  likelihood 
function  p(p|x)  should  be.  Precisely  defined,  the  likelihood  or  likelihood 
ratio  function  is  a  Radon-Nikodym  derivative  of  the  observation  probability 
measure  with  respect  to  another  probability  measure,  usually  of  an  observation 
noise  [27].  In  our  case,  the  other  probability  measure  will  be  the  clutter 
measure  vt. 

Let  Pi(F|x)  denote  the  conditional  probability  that  the  observation  pt 
lies  in  the  measurable  set  F  of  M(Y)  given  that  the  true  state  x(t)  *  x,  and 
let  P(j(F)  denote  the  probability  that  the  Poisson  measure  vt  lies  in  F.  We 
are  looking  for  the  Radon-Nikodym  derivative  of  P^  with  respect  to  Pq.  That 
is,  the  likelihood  function  p(p|x)  is  a  real-valued  function  of  the  set  M(Y) 
of  counting  measures  p  and  the  set  X  of  states  x  which  satisfies  the  relation 

Pl(F|x)  -  /  p(5|x)  dP0(O  .  (7-13) 

F 

Because  the  distribution  Pi  is  uniquely  determined  by  its  Laplace  transform 
[26],  proving  Eq.  7-13  Is  equivalent  to  proving 

/  exp(-CU))  dPxU)  -  /  exp(-?U))  PU|x)  dP0(O  (7-14) 

for  all  nonnegative,  measurable  $.  Written  in  terms  of  expectations,  Eq.  7-14 
is  equivalent  to 

E  {exp(-ut  ( <p)  )  J  |  x(  t) }  =  E{p(vt|x(t))  exp(-vt(  $)  )|(t)  }  .  (7-15) 
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We  can  rewrite  the  left  side  of  Eq.  7-15  using  Eqs.  7-4,  7-5  and  the  Laplace 
transform  relations  (Eqs.  7-11  and  7-12).  Note  also  that  the  random  measure 
vc  is  independent  of  x(t),  and  the  random  variables  z^(t)  and  y^(t)  are  con¬ 
ditionally  independent  of  each  other  and  vt  given  x(t).  Thus,  we  have 

n 

E{exp(-pt(4>)  )  |x(t)  }  =  n  [l-Gi  +  GiHi(exp(-4>) )]  exp(-X(l-exp(-,j>))  )  . 

k=l 


(7-16) 


In  Eq.  7-16  we  have  suppressed  the  dependence  on  x(t)  of  the  conditional 


detection  probabilities  G^  and  the  conditional  observation  probabilities  %. 


We  will  continue  to  do  this  to  simplify  our  notation. 


Our  objective  is  to  find  a  function  p(p|x)  such  that 


n 

n  [l-Gi  +  GiHi(exp(-<Ji))]  exp[-X(l-exp(-(j>))  ) 
k~l  (7-17) 

-  e{p( vt |x( t) )exp(-vt(  $) ) |x( t)  } 


i 

B 


We  will  do  this  in  a  series  of  increasingly  more  complicated  results,  beginning 
with  the  case  of  a  single  target  with  perfect  detection  (i.e.,  for  which  n  =  1 
and  G;l(x)  =*  1  for  all  x) .  In  this  simple  case,  the  basic  ideas  are  clearest. 
The  treatment  of  imperfect  detection  is  a  simple  corollary.  Next  follows  the 
case  of  n  targets  with  perfect  detection  and  finally  the  general  case. 


T 

**  *■ 


Proposition  1 

Assume  that  the  conditional  probability  Hi  has  a  density  ki  with  respect 
to  the  intensity  measure  X;  that  is,  assume  that 


Hi(Bjx)  =  /  ki(y|x)  dX(y) 
Then  for  n  =  1  and  Gi  =*  1,  p(p|x)  is  given  by 

p(y|x)  -  p(ki)  . 


(7-18) 


(7-19) 


■*  V.  h  .  rf*.  V*.  WL  ”■  a  3 . 1 


■r.y: 


-.  a  -'a  -V  ^  •  t  - '  *  s  ■ 
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Note  that  Eq.  7-19  means  that 


n 

p(u|x)  =  /  ki(y|x)  dp(y)  =  l  kl(y(j)|x)  (7-20) 

j=l 

if  the  counting  measure  p  is  given  by  the  sum  in  Eq.  7-1. 

The  assumption  about  the  measurement  density  in  this  proposition  is 
basically  that  the  conditional  observation  probability  distribution  Hi  is  ab¬ 
solutely  continuous  [24], [25]  with  respect  to  the  clutter  measure  X.  In  more 
physical  terms,  we  are  assuming  that  no  measurements  will  occur  where  there 
are  not  some  false  alarms.  In  many  cases,  the  measures  Hi  and  X  will  have 
densities  hi  and  n  with  respect  to  Lebesgue  measure  defined  on  the  observa¬ 
tion  space  Y  =■  Rm-  For  example,  this  is  the  case  if  one's  observations  have 
continuous  values  and  false  alarms  occur  uniformly  over  the  area  of  sensor 
coverage.  Then  ki  is  given  by  the  ratio 

hi<ylx> 

M(y|x) - - -  ,  (7-21) 

n(y) 

provided  that  hi(y|x)  =*  0  whenever  n(y)  *  0  so  that  the  ratio  is  well-defined 
(this  condition  is  the  substance  of  the  assumption  in  Proposition  1). 

Proof  of  Proposition 

Suppose  that  4>( y )  and  ki(y|x)  are  bounded  as  functions  of  y  for  each  x 
and  define 

$e(y)  =  <Ky)  -  e  ki  (y  |x(t)  )  (7-22) 

for  all  real  c.  From  Eq.  7-11  and  the  independence  of  x(t)  and  note  that 

exp(-X(l-exp(-!j>e)) )  =  E{exp(-vt(<fre))|x(t)  }  .  (7-23) 
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Under  the  boundedness  conditions  we  can  differentiate  Eq.  7-23  with  respect 
to  e.  Doing  so  at  e  =  0  gives 

X(ki  exp(-<j>))  exp(-X(  l-exp(-<J>) )  )  =  E{vt(ki  exp(-vt(  <f>)))  |  x(  t)  }  (7-24) 

which  gives  the  desired  result  upon  noting  x(k^  exp(-<(i))  =  (exp(~4>)  ).  The 

result  for  unbounded,  nonnegative  <ji  and  k^  follows  from  standard  application 
of  Legesgue's  dominated  and  monotone  convergence  theorems  [24], [25]. 

Corollary  1 

Assume  the  same  hypothesis  as  in  Proposition  1  but  do  not  assume  that 
G^(x)  =*  1.  Then  p(p|x)  is  given  by 

p(p|x)  -  1  -  Gi(x)  +  Gx(x)  p(k!)  (7-25) 


Proof  of  Corollary  1 

Note  that  from  Eq.  7-24  in  Proposition  1  we  have 

G]Hi(exp(-<t.))  exp(-X(l-exp(-<ji))  )  =  E{Givt(k1)exp(-vt(  <j>)  )  |x(t)  } 

(7-26) 


From  the  basic  formula  (Eq.  7-23)  we  have 

(1-Gi)  exp(-X(l-exp(-40) )  =  E  {( l-Gi)exp(-vt(  <J>) )  | x(  t)  }  .  (7-27) 


Adding  Eqs.  7-26  and  7-27  together  gives  the  desired  result  (Eq.  7-17)  for  the 
case  n=l. 

The  result  for  multiple  targets  is  more  complicated  and  requires  some 
preliminary  notation  and  definitions.  For  any  counting  measure  u  on  Y  define 
the  new  counting  measure  y(n)  on  Yn  as  follows.  For  each  bounded  real  func¬ 
tion  f  on  Yn  define  the  integral 
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p(n)(f)  =  l  on(ii,...,in)  f (y(ii) , • • • ,y(in) ) 


(7-28) 


where  the  sum  is  over  all  sequences  ii,...,in  of  integers  such  that  l<ij<n, 
and  an  =  0  if  ij  =  i^  for  j  *  k,  an  *  1  otherwise.  Note  that  ai(y(l))  =  1 

for  all  y(l).  The  sequence  y(k)  is  the  one  used  to  represent  p  in  Eq.  7-1. 

It  should  be  clear  that  the  definition  does  not  depend  on  the  order  of  the 

y(k) .  The  sum  in  Eq.  7-28  defines  the  integral  of  a  counting  measure  on  Y°. 
Note  that  the  expression  in  Eq.  7-28  can  be  interpreted  as  the  sum  over  all 
sequences  ix,...,in  such  that  no  ij  appears  twice  in  the  same  sequence.  Also 
note  that  p(l)  »  p. 

Finally,  if  fx,...,fn  are  function  of  Y,  define  the  product  by 


f1x...xfn(y1,...,yn)  -  fi(yi)...fn(yn) 


(7-29) 


We  are  now  ready  to  state  and  prove  the  results  for  multiple  targets,  first 
for  the  perfect  detection  case. 


Proposition  2 


Assume  that  the  conditional  probability  H*  has  a  density  ki  with  respect 
to  the  intensity  measure  \  for  each  i=l,...,n.  Assume  also  that  the  detection 
probability  Gi(x)  -  1  for  all  i  and  x.  Then  the  likelihood  function  p(p|x)  is 


given  by 


p(p|x)  =*  p(n)(k1x...xkn) 


(7-30) 


Note  that  Eq.  7-30  means  that 


p( P | x)  -  /  n  kj(y j | x)  dp(n>(yi . yn) 

j=l 

n 

=  l  OnUl.  ••  •  .in)  n  kj(y(ij)|x)  . 


(7-31) 
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Proof  of  Proposition  2 

We  need  to  prove  Eq.  7-17  Is  true  when  =  1  and  p(p|x)  is  given  by  Eq. 
7-30.  As  in  the  proof  of  Proposition  1,  it  suffices  to  prove  this  for  the 
case  when  <j>  and  are  bounded  functions  with  respect  to  y.  We  have  proven 
the  case  n  *  1  in  Proposition  1.  We  will  prove  the  general  case  by  induction 
on  n.  Thus,  assume  that  Proposition  2  is  true  for  n,  and  let  us  deduce  the 
result  for  n  +  1. 

The  result  for  n  states  that 
n 

n  [H^(exp(-$) )  ]  exp(-A(l-exp(-$) ) ) 
i=*l  (7-32; 

-  E{vt(n)(k1x..  .xkn)exp(-vt(<(>)  )|x(t)  } 

As  in  Proposition  1,  substitute  the  perturbation  $E 

krrt-l  (7-33] 

for  $  in  Eq.  7-32  and  take  the  derivative  with  respect  to  e.  This  gives  the 
following  result  when  e  »  0: 
n 

n  [Hi  exp(-<|>)  ]  exp[-A(l-exp(-$))  ) 
i-1 

“  E{vt(n)(k1x...xkn)  Vt(kiH-l)  exp[-vt(<f>)|x(t))} 
n  n 

-  I  A(kjkn+i  exp(-$) )  n  [Hi(exp(-<j>) )]  exp[-A(  l-exp(-<j>))  )  . 

j~l  i*j 

(7-34; 

Substituting  kjkfj+i  for  kj  in  the  result  (Eq.  7-32)  for  n  gives 

n 

*(k1kn+l  exp (-$))  n  [Hi  (exp(-$)  )  ]  exp[-A(  l-exp(-<t>))  ) 

=  E[vt(n)(k1x. .  .xk-ik^ix. . .  xkn)  exp(- vt(  4.)  )  I  x(  t)  }  . 


(7-35 
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Finally,  note  that  for  any  counting  measure  y  and  functions  kj  defined  on  Y 
defined  on  Y  we  have  that 


y(n+l)(k1x...xkrt+1)  =  y(n)(k1x...xkn)y(kn+1) 


-  I  p(n)(kix...xkjkn+ix...xkn) 

j=l 


(7-36) 


*  i. 

e 


Putting  Eqs.  7-35  and  7-36  in  Eq.  7-34  gives  the  desired  result  for  n  +  1. 

The  general  case  with  imperfect  detection  follows  Proposition  2  in  the 
same  way  that  Corollary  1  followed  from  Proposition  1.  The  notation,  however, 
is  more  complicated.  Adjoin  a  new  point  6  to  the  observation  set  Y  and  call 
the  extended  set  Y.  Extend  any  counting  measure  y  defined  on  Y  to  a  counting 
measure  jj  defined  on  JT  by  defining  _y(B)  if  8  is  not  in  B  and  y(B)  +  1  if  0  is 
in  B.  The  point  8  represents  a  fictitious  measurement  indicating  a  possible 
missed  detection  and  is  always  counted  by  _y.  Likewise,  extend  y(n)  to  _y(n) 
as  follows.  Define  Op  as  was  defined  before  except  that  ^  -  0  if  ij  -  i^ 
for  j  *  k  and  if  y(ij)  *  0.  Note  that  _y(n)  can  be  interpreted  similarly  to 
y(n).  The  sum  corresponding  to  Eq.  7-28  is  over  all  sequences  ij,...,in  such 
that  no  ij  appears  twice  in  the  same  sequence,  except  that  the  index  i  for 
which  y(i)  *  0  can  occur  any  number  of  times. 

It  is  also  necessary  to  extend  ki  to  a  function  k^  defined  on  Y  as 
follows : 

Jii(y | x)  =  Gi(x)  ki(y|x)  (7-37) 

for  y  in  Y,  and 

ki(Qjx)  =»  1  -  Gi(x)  .  (7-38) 
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Corollary  2 

Assume  the  same  hypothesis  as  in  Proposition  2  but  do  not  assume  that 
Gj(x)  =  1.  Then  p(p|x)  is  given  by 

p(y)x)  *  jj(n)(k2  x. . .  xkp)  .  (7-39) 

Proof  of  Corollary  2 

The  proof  follows  that  of  Proposition  2  with  minor  variations.  We  have 
proven  the  n=*l  case  in  Corollary  1.  To  prove  the  general  case  by  induction, 
assume  that  the  following  result  is  true  for  n  and  all  bounded  functions  4>,  kj 
and  detection  probability  functions  Gj. 
n 

n  [l-Gi  +  GiHi(exp(-$))]  exp (-X( 1-exp (- $) ) ) 
i-1  (7-40) 

-  E{vt(n)(k1x...xkn)  exp(-vt(4>))|x(t)  }  . 

Substitute  the  perturbation 

♦e  "  ♦  “  eGn+lknH  (7-41) 

for  <J>  in  Eq.  7-40.  The  derivative  at  e=0  gives 

n 

Gn+-lHn+l(exP(_<t>)]  n  [1_Gi  +  GiHi(exp(-<|>))]  exp[-X(l-exp(-<f.))  ) 
i=*l 

“  Ebit(n)(Jilx’-*xiSn)  Gn+1  v(^n+l)  exp(-vt(<j>)  ]  |x(t)  } 

(7-42) 

n  n 

_  I  Gn+lGjHj(kn+lexP(-<f>)  )  n  [l-Gi+GiHi(exp(-(j>)]] 

J-l  i*j 

exp(-X(l-exp(-^)) )  . 

Suppose  that  in  Eq.  7-40  we  substitute  the  function  f5j  =  Gn+iGjkn+^kj  for  kj 
and  1  for  Gj.  Then  we  obtain 
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E {_vt ^ n) (k_i  x . .  .x£jx. .  .xkp)  exp(-vt(  <J>)  )  |x(  t)  } 

n 

=  Gn+1GjHj(klH.1exp(-4i) )  II  [l-Gi  +  G^Hi  (exp(-$)  )  ]  exp(-A(l-exp(-<t,))  ) 

i*j 

( 

Multiply  Eq.  7-40  by  1  -  Gn+i  to  get 


(7-43 


(l“Gn+l)  11  [1_Gi  +  GlHi (exp( -$))  ]  exp(-A(l-exp(-<{>))  ) 

1=1  (7-44 

=  E(_ut(n)(kix...xkn)(l-GIt+i)  exp(-vt(<j»)  )|x(t)  }  . 


Adding  Eq.  7-44  to  Eq.  7-42  and  using  the  relation  7-43  gives 


n  [l-G*  +  G1Hi(exp(-4>))  J  exp (-A(  1-exp (-$))) 

1=1 

-  E{[vt(n)(k1x...xkn)vt(kxri.1) 
n 

“  I  vtCn)(k1x...x0jx...xkn)  ]  exp(-vt( $)  )  |x(t) } 

j-1 


(7-45 


(7-45 


Finally,  let  us  show  that  the  right  side  of  Eq.  7-45  is  what  we  want.  Corre¬ 
sponding  to  Eq.  7-35  is  a  recursive  expression  for  given  by 


n 

.  (il* • • • »*-n+l )  *  jSn(il » •  •  •  **n)  1  “  I  J§(i j»in+l ) 

L  i-1  J 


(7-46 


where  6(i, j)  =  1  only  If  l=j  and  y(i)  *  0  (j5  is  0  otherwise).  This  recursion 
implies  that  for  any  counting  measure  p  defined  on  Y  and  functions  Jy  defined 
on  the  extended  set  Y  we  have 

n 

p( n+1 ) (ki  x . . . xkn+i  )  =  _p(n)(ki  x. . .  xkn)  p(kp+i  )  -  £  jj(n) (k^  x. . .  xgj  x. . .  xkp) 

j=l 

(7-47 
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where  0j  is  defined  so  that  J3j(y)  =  j^j(y)^n+l(y^  for  y  in  Y  and  0j(  0)  =  0. 

If  jcj  and  k^+i  are  defined  as  in  Eqs .  7-37  and  7-38,  then  j3j  in  Eq.  3-37  is 
precisely  the  extension  of  0j  in  Eq.  7-43  to  Tf.  Thus,  substituting  Eq.  7-47 
with  y  »  vt  into  Eq.  7-45  gives  the  desired  result. 

7.4  ESTIMATION  ALGORITHM 

The  conditional  probability  distribution  of  the  state  variable  x(t), 
given  all  observations  up  to  time  period  t,  contains  all  the  statistical 
information  about  x(t)  that  exists  in  the  measurements.  Let  xc(x)  denote  the 
density  of  the  conditional  probability  distribution  with  respect  to  some  mea¬ 
sure  y  (y  might  be  Lebesgue  measure,  for  example).  Similarly,  let  *t+l,t(x) 
denote  the  conditional  probability  density  of  the  state  x(t+l)  given  all  ob¬ 
servations  up  to  the  time  t.  The  recursive  procedure  for  computing  xt+l(x) 
from  a  new  observation  and  the  predicted  density  Ht+l,t(x)  is  Bayes'  rule. 

In  terms  of  a  likelihood  function  p(ytlx)  for  the  measurement  model,  Bayes' 


rule  is 


*t+l(x) 


"t+l.tC*)  P( Mt lx) 

/  *t+l,t(0  P(MtU)  dY(5) 


(7-48) 


In  this  section  we  sketch  briefly  how  Bayes'  rule  applies  to  the  likelihood 
function  of  Eq.  7-39  in  the  previous  subsection  leads  directly  to  the  measure¬ 
ment  update  computation  found  in  hypothesis  tree  type  multiobject  tracking 
algorithms . 

Hypothesis  tree  type  estimation  algorithms  are  based  on  the  representa¬ 
tion  of  the  conditional  probability  density  of  the  state  as  a  weighted  sum  of 
simpler  densities,  namely. 
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Mx)  -  I  flj  ifj(x) 
j 


(7-49) 


where  each  itj(x)  is  a  product 


nj(x)  =  n  Xi^jCxi) 

i=l 


(7-50) 


The  density  1,i>j(xi)  represents  the  conditional  probability  of  the  ic^  target 
being  in  state  x^  given  that  hypothesis  j  is  true.  The  weight  0j  is  the 
probability  that  hypothesis  j  is  true.  For  example,  in  applications,  the 

densities  *i,j  and  xj  might  be  approximated  by  Gaussian  distributions  so  that 

Gq.  7-49  is  a  Gaussian  sum  distribution. 

Using  the  likelihood  function  of  subsection  7.3,  we  will  derive  the 
expression  (Gq.  7-49)  for  *t+l(x).  First  assume  that  x(t)  is  the  joint  vector 
of  individual  target  states  xj(t),  l<i<n.  For  simplicity  we  will  assume  that 
there  is  a  known  finite  number  n  of  targets  and  that  the  processes  x^(t)  are 
independent  Markov  chains.  Assume  that  the  measurement  Pt  *s  modeled  as  in 

subsection  7.2  and  that  the  measurement  densltites  hi(y|x)  satisfy 


hi(y|x)  -  h(yjxi)  , 


(7-51) 


and  the  detection  probabilities  Gi(x)  satisfy 


Gi(x)  =  G(xi) 


(7-52) 


Let  n(y)  denote  the  false  alarm  density,  and  as  in  subsection  7.3,  assume  that 
n(y)  is  not  0  unless  h(y|xi)  is  0  for  all  values  of  x^ .  In  this  case,  k^  is 


given  by 


G(*i)  h(y|xi) 
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for  y  in  the  observation  space  Y  (i.e.,  y  corresponding  to  a  detection),  and 
is  given  by 

k*(y|x)  =  k(y|xi)  =  1  -  GCxi)  (7-54) 

For  y  =•  8  (i.e.,  for  y  corresponding  to  a  missed  detection).  The  algorithm 
for  updating  Eq.  7-49  is  then  given  in  terms  of  a  prediction  step  and  a  mea¬ 
surement  step  as  follows. 

Prediction  Step 

For  each  hypothesis  index  j  and  each  target  index  i,  predict  ahead  the 
density  one  time  period  to  obtain  the  conditional  density  r;[jj+(x^)  of 

xi(t+l)  being  in  state  xj  given  measurements  up  to  time  period  t.  The  pre¬ 
dicted  density  is  given  by  the  weighted  sum  of  products 


"t+l.tCx)  »  l  8j  n  "i ,j+(*i) 
j  i-1 


(7-55) 


For  linear  Gaussian  target  models,  this  step  involves  just  the  prediction  of 
the  conditional  mean  and  covariance.  For  nonlinear  target  models,  this  step 
would  require  propagation  of  some  other  statistics  describing  the  conditional 
density  or  an  approximation  (such  as  in  an  extended  Kalman  filter).  Note  that 
the  prediction  step  does  not  change  the  hypotheses  j  or  the  probabilities  Bj 
of  the  hypotheses.  This  occurs  in  the  next  step  when  measurements  at  time 
period  t+1  are  processed. 

To  work  out  the  measurement  processing  step,  one  substitutes  the  expres¬ 
sion  7-55  for  *t+l,t(x)  and  the  likelihood  function  7-39  for  p(ut+llx)  into 
Bayes'  rule  (Eq.  7-48).  Let  us  carry  this  out  in  detail.  Suppose  that  the 
point  process  measurement  Mt+l  consists  of  the  observations  y(k),  l*k<m, 
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listed  in  some  arbitrary  order  (compare  Eq.  7-1).  Let  y(mfl)  denote  the 
fictitious  observation  y(nrt-l)  =  0.  Then  the  likelihood  function  p(Ut+llx) 
Eq.  7-39  is  given  by  the  summation 


P<ut+llx)  =  JJt+l^nHkix..  .kp) 

“  I  ^(il.-'-.in)  k(y(i1)|x1)...k(y(in) jxn) 


(7-56) 


where  the  sum  is  over  all  sequence  ij,...,in  of  integers  such  that  l<ij<m+l, 
and  Op  -  0  only  if  ij  ■  i^  for  j  *  k  and  ij  *  m+1,  else  op  *  1.  Note  that 
a  sequence  i},...,in  associates  target  j  with  one  of  the  m  observations  y(ij) 
or  with  the  fictitious  measurement  9  (to  denote  the  lack  of  real  observation 
for  that  target).  The  coefficient  oip  indicates  which  associations  are  valid, 
namely,  those  which  associate  each  real  measurement  with,  at  most,  one  target. 

Before  presenting  the  update  algorithm,  it  is  necessary  to  introduce  some 
further  notation.  If  measurement  k  is  associated  with  target  i  of  hypothesis 
j,  then  the  conditional  density  *i,j+  is  updated  to  xj,  j  by 

*i, j+(xi)G(xi)h(y(k) |xi) 

ni, j,k(xi)  =  - -  (7-57) 

/  iri,j+(5i)G(Ci)h(y(k)|^i)  dyUi) 

if  l<k<m,  and 

*i, j+(xi)[l~G(xi)] 

xi, j,nrt-l(xi)  “  - -  (7-58) 

J  xi,  j+(  5i)  [ l“G(  Ci)]  dy(0 

for  k  *  rnt-1  (processing  the  missed  detection).  Similarly,  the  likelihood  that 
the  measurement  k  is  associated  with  target  i  of  hypothesis  j  is  given  by 


°ij,k 


*i ,  j+( ^i)G( ?i)h (y(k) |  ) 

nfvrf’U'*  ) 


dY(Ci) 


(7-59) 
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If  l<k<m,  and  by 


°i ,  j ,m+l  =  /  n,j+(Ci)[l-G(5i)]  dY(5i)  (7-60) 


for  k  =»  nri-1  in  the  missed  detection  case. 


Using  this  notation,  substitute  Eqs.  7-56  and  7-57  to  obtain 


n 

*t+l,t(x)  P(Mtlx>  "  I  l  Sj  “n(il» •  *  *  »in)  n  ok.j.i  *k,j,i  (xk)  •  (7-61) 

j  i  k=l  x  k 


Thus,  Bayes'  rule  (Eq.  7-48)  gives 


n 

l  l  Bj  Op(ii,*..,in)  n  0k,j,i  ,xk,j,i,  <xk) 

4  4  1_^  I  fc,  Jt 

W»  -  - - - -  <-> 

1  I  Pj  °ji(^l » •  •  •  *^n)  n  °k,j,i, 
j  i  k-1  k 


The  update  step  of  the  estimation  algorithm  is  then  given  as  follows. 


Update  Step 

For  each  hypothesis  j,  each  target  index  i  and  each  measurement  index  k, 
update  the  density  to  obtain  the  updated  target  densities  xi,j,k(xi)  of 

x^(t+l)  being  in  state  xj  given  hypothesis  j  was  true  up  to  time  t  and  that 
measurement  k  was  associated  with  target  i  at  time  t+1 .  Similarly,  calculate 
the  likelihood  of  measurement  k  being  associated  with  target  i  given 

hypothesis  j  was  true  up  to  time  t.  For  each  hypothesis  j  at  time  t  and  for 
each  association  ii,...,in  of  measurements  with  targets,  generate  a  new 
hypothesis  j*  at  time  t+1  with  probability 


n 

On^l . !-n)  n  °k,j,i 

k=*l  x 


I  I  jjitil » •  •  •  »in)  n 

j  i  k-1  k 


(7-63) 
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and  with  target  conditional  densitities 

"k.j'Uk)  ■  *k,j,ik(xk)  •  (7-64) 

The  conditional  density  of  the  combined  target  states  at  time  t+1  is  then 

n 

*t+l(x)  “I  Gj*  n  .  (7-65) 

j’  i*l 

Note  the  tree  structure  of  the  hypotheses.  At  every  update  step  each 
hypothesis  j  generates  many  new  disjoint  hypotheses  j'.  Over  time,  the 
hypotheses  can  be  arranged  in  a  tree  with  the  nodes  at  one  level  of  the  tree 
representing  hypotheses  at  one  period  of  time .  Branching  of  nodes  at  one 
level  of  the  tree  to  nodes  at  the  next  level  of  the  tree  represents  the  gen¬ 
eration  of  new  hypotheses.  Note  also  that  the  number  of  nodes  in  this  tree 
increases  roughly  as  nmC  if  m  is  the  number  of  observations  per  time  period, 
t  is  the  number  of  time  periods,  and  n  is  the  (fixed)  number  of  targets. 
Several  approximations  have  been  developed  to  limit  this  geometric  growth  of 
hypotheses  (see  [1],[2])- 

7.5  APPLICATION  TO  PERFORMANCE  ANALYSIS 

In  this  section  we  discuss  briefly  how  the  result  of  subsection  7.3  can 
be  used  to  obtain  Cramer-Rao  type  lower  bounds  of  mean  square  terror  in  esti¬ 
mating  x(t)  given  a  measurement  model  of  the  form  described  in  Eq.  7-4  and 
7-5.  From  Section  4,  note  that  the  Cramer-Rao  bound  depends  on  the  Fisher 
information  defined  by 

!1  3pT  3p  ) 

-  —  }  •  (7-66) 

p(ut  |x(t)  )2  3x  3x  ) 
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The  derivatives  are  computed  from 

3p  n  f  *jkj  N 

—  (yjx)=  \  _p(n)  kj^x...x  —  x.^xj^  •  (7-67) 

3x  j=l  V.  3x  ) 

It  is  assumed  that  the  measurement  density  k^(y|x)  and  the  detection  probabil¬ 
ity  G^(x)  are  both  differentiable  with  respect  to  x.  The  expressions  needed 
to  compute  I(t)  reduce  to 

p(y|x)  =»  l-Gi(x)  +  Gi(x)  p(ki(y  |x)  )  (7-68) 

and 

3p  3G  3G  v  ( 3k  ^ 

—  (y|x)  -  -  —  l(x)  +  —  l(x)p(k1(y|x)  )  +  Gi(x)p  —l(y|x)  (7-69) 

3x  3x  3x  \3x  J 

for  the  case  of  one  target  (n»l). 

Actually  computing  I(t)  requires  taking  the  expectation  in  Eq.  7-66  with 
respect  to  the  random  variable  x(t)  and  the  random  measure  yt.  In  general, 
this  expectation  is  not  expressible  in  closed  form  and  it  is  necessary  to 
resort  to  analytical  (e.g-,  perturbation  expansion)  or  numerical  approximation 
(e.g.,  perturbation  expansion)  or  numerical  approximation  (e.g.,  Monte-Carlo 
simulation).  Note  that  once  I(t)  has  been  obtained,  it  is  easy  to  find  the 
Cramer-Rao  lower  bound  on  the  error  covariance  for  estimating  a  Gaussian  dis¬ 
tributed  state.  Suppose  that  x(t)  satisfies  the  linear  Gaussian  equation 

x(t+l)  ■  A(t)x(t)  +  B(t)w(t)  (7-70) 

where  w(t)  is  Gaussian  white  noise  with  covariance  Q(t)  and  x(l)  has  prior 
covariance  Iq*  Initialize  Pq  and  So  to  Eq,  and  compute  Pt  and  St  using  the 
recursion  defined  by 
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pt+l  -  [st_1  +  I(t+1)  ]_1  (7-71) 

for  t>0  and 

St  -  A(t)PtA(t)T  +  B(t)Q(t)B(t)T  .  (7-72) 

for  t>l.  If  x(t)  is  the  minimum  mean  square  estimate  of  x(t)  based  on  the 
measurements  ps  for  l<s<t,  then  Pt  is  the  Cramer-Rao  lower  bound  on  the  error 
covariance 

Pt  <  Ej[x(t)-i(t)][x(t)-x(t)]T}  .  (7-73) 

7.6  CONCLUDING  REMARKS 

In  this  section  we  have  introduced  a  new  formulation  of  multiobject 
tracking  as  an  estimation  problem  with  random  point  process  measurements  and 
derived  the  corresponding  likelihood  ratio  formula  using  Laplace  transform 
methods.  This  new  theoretical  approach  provides  a  compact  formulation  and 
powerful  technique  for  analyzing  multiobject  tracking  problems.  It  allows 
easy  derivation  of  optimal  tracking  algorithms  and  may  prove  useful  in  ana¬ 
lyzing  the  optimal  performance  of  such  algorithms.  The  approach  also  extends 
easily  to  other  kinds  of  multiobject  tracking  problems  than  the  one  described 
in  subsections  7.3  and  7.4,  such  as  measurement  models  in  which  several  re¬ 
turns  are  received  from  the  same  target  (e.g.,  multipath  phenomena).  Although 
much  work  remains  to  be  done,  we  believe  that  the  approach  presented  in  this 
section  provides  an  effective  framework  for  developing  new  analytical  tools 
to  understand  better  the  quantitative  performance  of  multiobject  tracking 


algorithms. 
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SECTION  8 

NUMERICAL  EXAMPLES  AND  CONCLUSIONS 

8.1  NUMERICAL  COMPARISON  OF  PERFORMANCE  METHODS 

To  understand  better  the  relative  accuracy  of  the  different  performance 
methods  presented  in  Section  4  through  6,  we  considered  the  following  simple 
Type  I  hybrid  state  estimation  problem: 


x(t+l) 

-  X<t)  +  w(t) 

(8-1) 

y(t) 

-  x(t)  +  q(t)vx(t) 

(8-2) 

+  (l-q(t)  )v0(t)  . 


This  is  a  dynamic  version  of  the  simple  static  problem  analyzed  in  Section  3. 

We  assume  that  the  variance  oi^  of  v^(t)  Is  much  larger  than  the  variance  oq 2 
of  vq(0-  That  is,  q(t)  =  1  indicates  a  bad  measurement  and  q(t)  =  0  indicates 
a  good  measurement • 

We  describe  two  cases  whose  parameters  are  shown  in  Table  8-1.  The  two 
cases  differ  only  by  the  initial  variance  of  x(0).  The  initial  variance  is 
small  in  Case  1  to  let  us  study  steady-state  estimation  performance  —  perfor¬ 
mance  achieved  once  the  initial  uncertainty  in  the  state  estimate  has  been 
reduced.  The  large  initial  variance  in  Case  2  lets  us  study  acquisition 
performance  —  the  details  of  how  a  very  large  initial  uncertainty  is  reduced 
until  steady-state  performance  is  achieved. 
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TABLE  8-1.  STATISTICAL  PARAMETERS  FOR  NUMERICAL  EXAMPLES 


CASE  1 

CASE  2 

Var  (x(0)) 

4.0 

104 

Var  (w) 

4.0 

4.0 

Prob  (q=l ) 

.1 

.1 

Var  (vjJ 

104 

104 

Var  (v0) 

4.0 

4.0 

Figures  8-1  and  8-2  show  performance  estimates  for  the  two  cases.  In 
each  figure  we  have  plotted  four  types  of  performance:  Cramer-Rao  lower  bound 
(CRB),  an  upper  bound  given  by  the  best  linear  estimator  (LIN),  simulated  per¬ 
formance  of  a  simple  gating  algorithm  (GATE) ,  and  simulated  performance  of  the 
optimal  algorithm  (OPT).  We  also  computed  the  rate  distortion  lower  bound  of 
Section  5  for  this  problem,  but  the  performance  prediction  was  several  orders 
of  magnitude  worse  than  the  other  methods,  so  we  omitted  the  rate  distortion 
results. 

Figure  8-1  shows  steady-state  performance.  The  best  linear  estimator's 
performance  increases  since  it  assumes  measurement  noise  is  Gaussian  with 
variance  equal  to  the  average  variance  of  vi(t)  and  vg(t), 

e  +  (v-e)oo2  =  (,1)104  +  (.9)4 

=  103  . 

As  t  +  °»,  the  best  linear  estimate  will  achieve  steady-state  performance,  but 
only  for  t  much  greater  than  the  seven  time  steps  considered  in  Fig.  8-1. 
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Figure  8-1.  Case  1:  Steady-State  Performance 
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TIME 

Figure  8-2.  Case  2:  Acquisition  Performance 
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The  Cramer-Rao  lower  bound,  the  simulated  gate  algorithm,  and  the  simu¬ 
lated  optimal  algorithm  all  track  fairly  well  in  Fig.  8-1.  Note  that  one 
expects  the  gate  algorithm  to  be  near-optimal  in  steady-state:  that  is,  it 
performs  well  once  a  good  state  estimate  has  been  achieved. 

Figure  8-2  shows  performance  when  the  initial  state  variance  is  the  same 
magnitude  as  bad  measurements.  Expressed  in  tracking  terms,  an  estimator  must 
first  acquire  the  target:  that  is,  reduce  the  initially  large  error  variance 
to  steady-state  performance.  As  shown  in  Fig.  8-2,  the  optimal  estimator 
acquires  the  state  after  about  three  measurements.  On  the  other  hand,  the 
gate  algorithm  never  acquires  the  target  and  the  best  linear  estimator  does 
so  only  very  slowly.  The  Cramer-Rao  lower  bound  optimistically  predicts  a 
quicker  acquisition  (after  1  measurement),  but  predicts  steady-state  perfor¬ 
mance  accurately.  Note  that  the  poor  performance  of  the  gate  algorithm  is  to 
be  expected.  Gate  algorithms  are  designed  based  on  steady-state  assumptions 
and  usually  require  an  auxiliary  initialization  algorithm  to  acquire  the 
steady-state  performance. 

As  a  final  numerical  example,  Table  8-2  shows  results  comparing  the 
Monte-Carlo  method  of  Section  7  with  the  simulated  optimal  algorithm  for  Case 
2  above  over  the  first  five  measurement  updates.  We  have  considered  three 
different  choices  of  the  Monte-Carlo  parameters  n  and  m  (see  Section  7).  Note 
that  the  Monte-Carlo  performance  method  predicts  the  acquisition  phenomena  but 


tends  to  be  pessimistic  about  steady-state  performance. 
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TABLE  8-2.  MONTE-CARLO  PERFORMANCE  METHOD 


8.2  CONCLUDING  REMARKS 

Analyzing  the  optimal  performance  achievable  In  multiobject  tracking  pre¬ 
sents  two  difficulties:  the  estimation  problems  are  Inherently  non-Gausslan 
and  the  dimensions  are  very  large.  In  Sections  3  through  6  we  addressed  the 
non-Gausslan  problem  and  In  Section  7  we  addressed  the  dimensionality  Issue. 

Of  the  different  methods  we  Investigated,  three  appear  promising  enough  to 
warrant  further  research:  the  Cramer-Rao  method  of  Section  4,  the  represen¬ 
tation  theorem  method  of  Section  6,  and  the  random  point  process  approach  of 
Section  7.  The  Cramer-Rao  method  of  Section  4  has  the  attractive  computational 
property  that  Its  complexity  only  decreases  linearly  with  time.  Furthermore, 
numerical  examples  suggest  the  lower  bound  is  an  accurate  approximation  of 
steady-state  performance  in  some  cases.  We  need  to  understand  now  precisely 
what  these  cases  are,  and  prove  the  validity  of  the  Cramer-Rao  bound  as  a 
steady-state  approximation. 


ALPHATECH,  INC 


The  Monte-Carlo  approach  of  Section  6  is  promising  but  requires  further 
work  before  we  can  decide  conclusively  whether  or  not  it  is  practically  useful. 
The  advantage  of  this  approach  is  that  it  can  give  an  accurate  estimate  of 
optimal  performance,  whether  in  acquisition  or  steady-state  phase  —  if  one 
computes  long  enough.  The  disadvantage  is  that  this  may  require  impractically 
long  computation.  The  computational  aspects  of  the  problem  require  further 
investigation  and  we  indicated  some  promising  directions  in  Section  6. 

In  Section  7  we  presented  a  new  theoretical  approach  to  analyzing  the 
type  of  hybrid  state  problem  arising  in  multiobject  tracking.  We  also  sketched 
how  this  approach  might  be  used  with  Cramer-Rao  methods  to  obtain  performance 
bounds.  We  believe  the  random  point  process  formulation  and  techniques  of 
Section  7  will  prove  useful  in  analyzing  multiobject  tracking  problems  by 
other  methods  as  well.  This  approach  seems  especially  fruitful  for  seeking 
analytical  results  as  the  formulation  allows  one  to  represent  complex  multi¬ 
object  tracking  problems  very  compactly. 
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APPENDIX  A 


SECOND-ORDER  EXPANSION  OF  MINIMUM  MEAN  SQUARE  ERROR 


As  shown  in  subsection  3.2.2  the  minimum  mean  square  error  V  is  given  by 


o2  o02 


°2  +  °o2 


o^[oi2  -  oo^]  c i^[oi2  -  oo^  ]2 

[c2+0()2]2  '  [ a2  +  oq2  ]2  [a2  +  0l2]2 


(A-I) 


where 


-  52  4>lU)2 

J  -  /  -  , 

-»  ( 1— e)  4>o( C)  +  e  4>l(0 


(A-2) 


4»iU)  -  (2*  Qi)-1/2  exp 1/2  52  Qi"1 }  , 


(A-3) 


Qi  -  o2  +  0i2  . 


(A- A) 


We  wish  to  approximate  J  in  the  case  that  a  "  oj  »  oq  and  e  is  small. 


Reduction  of  Integral 


By  making  the  substitution 


we  find  that 


J  *  a  I 


(A-5) 


where  I  is  the  Integral 


‘•1>J 
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i 


j"  exp{-C2} 

-<*>  1  +  p  exp  {-Xc^  } 


4  e-1 


(A-6) 


(A-7) 


(A-8) 


Using  Che  substitution 


u  ■  X^2  -  £np 


gives  us 


I  ■  -  yl/2  |£npjl/2  p~Y  .  k 


(A-9) 


where  the  integral  K  is 


K 


.  f  1  +  -a-  11/2  o_UT 

l  frip  ; _ 

-  £np  1  +  e“u 


and 


Y  =  V1 


(A-10) 


The  integral  K  can  be  written  as  the  sum 


K  =■  Ki  +  K2 


(A-ll) 


>  W.W. 


:  >iw 


.  .-V 


.-AS  .  .  ,  . 

.1^  --- 


,  - 1 
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and  K2  is  the  integral 


-inp  1  +  e-u 


u  W2  p-l 


The  second  integral  K2  can  also  be  expressed  as 


1  (1-z)1^2  *  e2*9 

K2  *  inp  •  /  — — -  d2 

0  1  +  e~2ZnP 


where 


0  »  (l-Y)tnp  . 


(A-12 


In  the  case  of  interest  to  us,  y  "  1,  in p  »  1,  but  9  is  of  order  1.  We  will 


now  approximate  Kj  and  K2  under  these  conditions. 


Approximation  of  K]  and  K; 


Write 


a  =»  Jtnp  , 


(A-13 


so  that  K^and  K2  become 


u  W2  .-1 


Kl  -  / 


00  1  -t-  e“u 


0  (l-z)l/2  e02 

K2  =  a  •  /  - *  dz 


l  +  e“z 
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where  we  assume  a  >  0,  y  *>  1,  and  0  ~  1.  We  approximate  as 


K1  =  K11  +  *12 


where 


"  e"Tu 

Kn  =  8(y)  =  /  -  du 

0  1  +  e_u 


is  the  incomplete  beta  function  and 


*12  = 


oo 

/ 

0 


1  +  e_u 


1 


du 


The  second  term  Kj2  can  be  bounded  by 


0  <  K12  <  —  /  ue-Yu  du 
2a  0 


1 

2a  y2 


Note  that 


BC1)  -  £n2  , 


e'd) 


“IT2 

12 


and 


0"(1)  <  /  u2e_Tu  du 
0 


2 


TT2  1 

6(y)  =  £n2  -  —  (y-1)  +  —  A(y) 


(1-Y)2 


Thus ,  we  have 
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where 


2 

&( y)  <  —  if  y  <  1 
Y3 

A(y)  <2  if  y  >  1  • 


Thus,  we  approximate  by 

it2 

Ki  3  £n2  -  —  (y~1) 


where  the  magnitude  of  the  error  is  bounded  by 

-  +  (1-Y)2  max{Y~3,  l}  • 

2a  Y2 

For  K2  we  use  the  crude  bounds 

1  i  2 

-  <  a  <  -  e® 

3  3 

if  8  >  0,  or 

—  e®  <  a_1K?  <  —  . 

3  3 

if  0  <  0. 

Upper  and  Lower  Bounds 

Let  us  assume  y  >  1;  this  is  equivalent  to  assuming 

oi2  -  02  <;  2oo2 

Furthermore,  assume  that  Jlnp  >  0.  This  is  equivalent  to 


(A-14) 


(A-15) 


(A-16) 


(A-17) 


(A-18) 
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o2  +  o02 


0^[oi2  -  OQ^  ]  0^[oi2  -  OQ^  ] 


[o2  +  0Q2 j  [a2  +  a02]2  [a2  +  o^2]2 


Asymptotic  Behavior  as  e  *  0 

Assume  that  y  >  1  is  maintained  as  before  and  let  e  ♦  0.  Then  note  that 


1  (l-z)l/2  top  e"z62nP  dz 
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Furthermore, 


Note  that  as  e  0 


lim  Ki  -  8(y) 
e+0 


-er  • 


lap  ~  In  £  , 


a  ~  e-1  • 


It  follows  that  as  e  +  0 


J  —  c  •  j  Jin 


where 


L2i  y3/2  f  —  ^ 

nr  \Qo  J 


Ql  V  l/2y 


In  this  case,  J  does  approach  0  as  e  ♦  0  but  not  very  rapidly  since  fl  is 
assumed  small. 
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